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INTERIOR  FLUID  DYNAMICS 
OF  LIQUID-FILLED  PROJECTILES 


1.  Research  Objectives 

The  project  "Interior  Fiuid  Dynamics  of  Liquid-Filled  Projectiles"  under  Contract 
DAAA15-85-K-0012  of  the  Department  of  the  Army,  AMCCOM,  was  originally  planned  as  a 
three-year  effort  with  the  working  period  07/01/85  -  06/30/88  and  Thorwald  Herbert  as  principal 
investigator.  The  working  period  was  extended  at  no  cost  until  12/31/88. 

A  detailed  description  of  the  research  objectives  has  been  given  in  section  2.2  of  the  original 
proposal.  In  summary,  the  research  aimed  at  understanding  the  origin  of  the  flight  instability  of 
projectiles  with  liquid  payloads  and  at  estimating  and  accurately  calculating  the  liquid-induced 
moments  to  obtain  the  data  necessary  for  the  design  of  future  projectiles.  This  effort  was  directly 
related  to  the  experimental  studies  conducted  at  CRDEC  and  to  numerical  studies  performed  at 
BRL.  In  detail,  the  work  under  this  contract  was  planned  to  extend  our  earlier  analysis  in  three 
directions: 

(1)  Fully  exploit  the  potential  of  the  linearized  solution  with  regard  to  yaw  moment,  pitch 
moment,  and  temperature  effects  on  the  viscosity.  Extend  the  solution  to  unsteady  (e.g. 
spin-up)  situations. 

(2)  Perform  a  perturbation  analysis  of  the  nonlinear  effects  that  produce  radial  and  azimuthal 
velocity  components  and  modify  the  axial  component.  Study  the  effect  of  the  cylinder  end 
walls  on  the  flow  field  and  the  pressure.  Investigate  the  onset  and  nature  of  cellular  motions 
at  higher  Reynolds  numbers. 

(3)  Develop  a  highly  efficient  computer  code  for  calculating  the  steady  flow  field  and  .Iquid 
moments  on  the  basis  of  the  Navicr-Stokes  equations.  Implement  spectral  approximations  in 
all  space  variables.  Extend  the  analysis  to  unsteady  problems  (spin-up,  spin-down,  yaw-angle 
growth).  Incorporate  the  flow  computations  into  an  existing  six-dcgrec-of-frccdom  code. 

As  will  be  discussed  in  the  next  section,  most  of  these  goals  have  been  reached.  Various 
other  results  have  been  achieved  that  were  not  anticipated  at  initiation  of  this  contract.  The  work 
on  unsteady  aspects  of  the  problem  was  initiated  but  not  finished  within  the  period  of  this  con¬ 
tract. 

Owing  to  Th.  Herbert’s  accepting  a  new  position  at  The  Ohio  State  University  (OSU)  begin 
ning  with  the  academic  year  1987/88  on  October  1,  1987,  only  a  part  of  the  program  was  con¬ 
ducted  at  the  Virginia  Polytechnic  Institute  and  State  University  (VPI).  A  part  of  the  remaining 
funds  was  made  available  through  a  subcontract  to  The  Ohio  State  University  to  continue  the 
research.  The  research  equipment  acquired  within  this  contract  remained  at  VPI.  The  lack  of  this 
equipment  required  time-consuming  modifications  of  our  research  codes  and  caused  delays  in  the 
research  program. 


2.  Research  Achievements 


After  initiation  of  this  contract,  work  was  simultaneously  conducted  on  all  three  topics. 
Rapid  progress  was  made  until  summer  1987  in  close  cooperation  of  the  principal  investigator 
with  Rihua  Li  who  worked  as  a  postdoctoral  associate.  At  this  time,  Dr.  Li  accepted  a  similar 
position  at  the  University  of  Arizona  since  the  funds  available  at  OSU  were  insufficient  for  his 
remaining  in  the  program.  At  OSU,  M.  Sclmi  continued  Dr.  Li’s  work  after  some  training  period. 

Most  of  the  results  of  our  research  have  been  reported  in  publications  (see  section  4)  and 
more  detailed  papers  arc  in  preparation.  Therefore,  we  can  restrict  this  report  to  a  brief  overview 
of  the  achievements. 

Our  earlier  analytical  studies  for  cylinders  of  large  aspect  ratio  (Herbert  1986)  were  extended 
to  the  nonlinear  problem  by  use  of  a  perturbation  method.  The  small  expansion  parameter  is 
e  -  xsinO,  where  t  =  Q  the  nutation  rate,  w  the  spin  rate,  and  0  the  angle  between  the  two 
axes  of  rotation.  The  solution  up  to  second  order  is  obtained  in  closed  form,  while  the  third-order 
contributions  arc  determined  numerically.  Higher-order  terms  arc  small  such  that  good  estimates 
of  all  liquid  moments  can  be  obtained  from  the  closed-form  approximation.  The  effect  of  finite 
aspect  ratios  has  been  evaluated  by  comparison  with  results  of  the  spectral  code  (sec  below).  This 
evaluation  and  the  comparison  of  (light  simulations  based  on  analytical  and  numerical  results 
show's  that  the  estimates  arc  conservative  and  therefore  useful  for  preliminary  design  studies. 
Recently,  the  analysis  has  been  extended  to  obtain  results  for  partially  filled  cylinders  with  an 
axisymmctric  void,  for  cylinders  with  a  central  rod,  and  for  the  case  of  two  immiscible  fluids  of 
different  viscosities  and  densities  separated  by  a  cylindrical  interface. 

Various  studies  have  been  performed  to  evaluate  existing  Navicr-Stokcs  solvers  (Vaughn  ct 
al.  1985a,  Strikwcrda  &  Nagel  1985).  These  efforts  were  kindly  supported  by  M.  Nusca,  BRL, 
who  provided  data  for  comparison  and  the  latest  version  of  the  Sandia  code.  The  previous  codes 
are  primarily  based  on  finite-difference  approximations  with  relatively  coarse  grid  and  were 
designed  without  insight  into  the  nature  of  the  fluid  motion.  As  a  result,  these  codes  provide  solu¬ 
tions  only  for  relatively  small  Reynolds  numbers  and  have  to  compromise  between  numerical 
approximation  and  computational  expense.  Guided  by  our  analytical  work,  we  have  directed  our 
efforts  towaid  tailoring  the  Navicr-Stokcs  code  to  the  nature  of  the  fluid  motion  and  obtaining 
more  accurate  moments  from  a  given  numerical  approximation  for  the  flow  fields. 

Previous  numerical  w'ork  determined  the  moments  from  surface  integrals  over  normal 
stresses  and  shear  stresses  which  involve  the  pressure  and  velocity  gradients  at  the  boundary. 
However,  boundary  quantities  arc  obtained  with  numerical  errors  typically  much  larger  than  the 
velocity  in  the  interior.  We  have  developed  a  method  to  determine  the  moments  by  integration 
over  die  volume  of  the  cylinder.  This  integration  requires  only  the  knowledge  of  the  fundamental 
Fourier  component  of  the  axial  velocity  and  the  mean  (streaming)  component  of  the  azimuthal 
velocity.  The  pressure  and  velocity  gradients  arc  eliminated  from  this  formulation.  In  comparison 
with  tlic  traditional  surface  approach,  die  smoothing  by  integration  in  three  space  directions  and 
the  use  of  well-approximated  data  provides  moments  of  superior  accuracy  with  the  same  original 
data,  i.e.  with  the  same  amount  of  computer  time.  The  method  is  generic  in  ,he  sense  that  it  is 
v ai id  not  only  for  cylinders  but  as  well  for  arbitrary  closed  containers.  Moi.  er,  the  method  is 
tmt  only  use! ul  for  purely  numerical  work  but  especially  for  perturbation  approaches  that  provide 
the  fundamental  Fourier  component  of  the  axial  velocity  at  first  order  and  the  mean  component  of 
the  azimuthal  velocity  at  second  order. 


The  design  of  the  spectral  code  for  solving  the  Navicr-Stokcs  equations  was  guided  by 
analytical  studies,  the  anticipated  use  for  flight  simulations,  and  the  desire  to  keep  the  code  open 
tor  extension  to  unsteady  problems.  The  latter  demand  prevented  the  use  of  artificial  (non¬ 
physical)  time  dependence  to  obtain  the  steady  solution.  Since  the  structure  of  the  equations 
clearly  exhibits  the  dominance  of  the  linear  terms  at  the  usually  small  values  of  e,  a  good  approxi¬ 
mation  for  the  solution  can  be  obtained  by  solving  the  linear  system.  This  solution  is  improved 
by  iteration  to  account  for  nonlinearity.  Besides  the  linearization  in  e,  a  second  type  of  lineariza¬ 
tion  can  be  made  in  the  velocity  deviation  from  rigid-body  motion.  A  third  type  of  linearization  is 
about  a  known  solution  at  neighboring  parameters.  While  the  second  type  is  used  in  absence  of 
prior  knowledge,  the  third  type  is  especially  suited  for  continuation  of  the  solution  and  moment 
calculations  in  flight  simulations  which  arc  characterized  by  a  slow  variation  of  the  parameters 
along  the  trajectory. 

The  utility  of  the  spectral  code  for  fundamental  studies  of  the  fluid  motion  (Herbert  &  Li 
19X7)  and  its  efficiency  for  flight  simulations  (Herbert  1988)  has  been  verified.  The  flight  simula¬ 
tions  were  based  on  a  slightly  modified  version  of  the  code  developed  by  Vaughn  ct  al,  (1985b). 
More  recent  work  was  dedicated  to  the  range  of  high  Reynolds  numbers  and  to  the  viscous  effects 
on  the  resonance  with  inertial  waves.  For  sufficiently  small  e  the  spectral  solutions  of  the  non¬ 
linear  problem  agree  well  with  the  results  of  Hall  ct  al.  (1987)  obtained  by  spatial  eigenfunction 
expansions  of  the  linearized  problem. 

The  eigenfunction  expansions  used  by  Hall,  Sedney  &  Gerber  (1987)  require  numerical 
determination  of  eigenvalues  and  eigenfunctions  from  an  eigenvalue  problem  for  a  system  of  ordi¬ 
nary  differential  equations.  We  have  developed  an  equivalent  but  much  simpler  set  of  equations 
that  permits  closed-form  solutions  for  the  spatial  eigenfunctions  (Li  &  Herbert  1988)  and  thus 
avoids  the  high  computational  expense  for  generating  numerical  solutions.  The  eigenvalues  are 
obtained  as  solutions  of  a  transcendental  equation.  Our  formulation  of  the  problem  clearly  reveals 
the  structure  of  the  eigenvalue  spectrum  and  explains  the  empirical  results  of  Hall  ct  al.  on  the 
grouping  of  the  eigenfunctions.  The  eigenvalues  have  been  generated  for  Reynolds  numbers  as 
high  as  106  The  numerical  problems  of  using  modified  Bessel  functions  of  large  complex  argu¬ 
ments  at  these  high  Reynolds  numbers  arc  not  yet  fully  overcome. 

Besides  solving  the  new  formulation  in  terms  of  eigenfunction  expansions,  we  apply  spectral 
expansions  in  Chebyshev  series.  The  otherwise  two-step  solution  procedure  (generation  of  func¬ 
tions,  solving  for  the  expansion  coefficients)  reduces  in  this  case  to  a  single  step.  Although  this 
code  has  not  yet  matured  for  routine  applications,  experience  in  the  range  up  to  Reynolds  numbers 
of  the  order  of  103  is  very  encouraging. 

Various  studies  have  been  performed  to  adapt  the  existing  analytical  and  numerical  methods 
io  ihe  unsteady  problem.  These  studies  concentrate  on  more  efficient  solution  of  linear  and 
weakly  nonlinear  algebraic  systems  with  iterative  methods. 

With  the  work  under  this  contract,  we  have  deepened  flic  understanding  of  the  fluid  motion 
in  liquid-filled  cylinders.  The  analytical  work  provides  estimates  for  the  moments  in  various 
cylindrical  configurations,  guidance  for  the  design  of  numerical  methods,  and  an  improved  basis 
for  calculation  of  the  moments.  Spectral  codes  for  solving  the  linear  and  nonlinear  problem  have 
been  developed  and  verified  to  efficiently  calculate  the  liquid  moments  up  to  Reynolds  numbers  of 
flic  order  ol  210\  The  spectral  code  for  the  linear  problem  will  lie  further  developed  to  overlap 
with  the  boundary-layer  methods  in  the  range  of  high  Reynolds  numbers. 


3.  Personnel 


During  the  working  period,  the  following  personnel  were  partly  supported  under  Contract 
DAAA15-85-K-0012: 

Thorwald  Herbert,  Professor,  Principal  Investigator 

Rihua  Li,  Postdoctoral  Associate 

Rclja  Zivojnovic,  Graduate  Student  (M.S.  level) 

Stephen  D.  Greco.  Graduate  Student  (Ph.D.  level) 

Mohamed  Sclmi,  Graduate  Student  (Ph.D.  level) 

Charlotte  R.  Hawley,  Research  Specialist 
Vinccl  Menta,  Undergraduate  Student 
David  Pierpont,  Undergraduate  Student 

Rihua  Li  developed  the  volume  approach  for  the  calculation  of  the  liquid  moments  and 
implemented  this  method  into  the  numerical  work.  L.i  analyzed  the  problems  associated  with 
comers  of  the  computational  domain  for  spectral  codes  and  contributed  both  to  the  analytical  and 
numerical  aspects  of  the  research  program.  In  1987,  Li  developed  the  simplified  formulation  of 
the  linear  problem  and  the  closed-form  solution  for  the  eigenfunctions.  Dr.  Li  is  presently  Post¬ 
doctoral  Associate  in  the  Department  of  Aerospace  and  Mechanical  Engineering  at  the  University 
of  Arizona. 

R.  Zivojnovic  and  S.  Greco  were  involved  in  developing  the  perturbation  approach  and  solv¬ 
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Abstract 

In  the  framework  of  a  feasibility  study,  we  have  designed  a  small  model  test  fixture  for  visualization 
of  the  flow  in  a  spinning  and  nutating  cylinder.  We  describe  the  apparatus  and  the  visualization  tech¬ 
nique,  and  report  some  results.  As  the  Reynolds  number  increases,  we  observe  an  axially  almost  uniform 
flow  that  turns  at  the  ends,  the  development  of  two  elongated  cells  in  the  plane  of  the  spin  and  nutation 
axis,  the  formation  of  additional  laminar  cells,  and  ultimately  unsteady  and  turbulent  flow  with  a  super¬ 
posed  large-scale  cellular  motion. 


1.  Introduction 

It  is  well-known  that  spin-stabilized  shells  carrying  liquid  payloads  can  suffer  dynamical  instability. 
For  cylindrical  cavities  and  low  viscosity  of  the  liquid,  the  instability  due  to  basically  inviscid  inertia! 
waves  is  rather  well  understood1.  The  instability  of  certain  shells  like  the  XM761,  however,  is  dis¬ 
tinguished  in  character  by  the  rapid  loss  in  spin  rate.  Experiments2  and  subsequent  field  tests3  establish 
that  this  flight  instability  is  most  pronounced  for  liquid  fills  of  very  high  viscosity. 

Theoretical  analysis  of  a  simple  model  of  the  internal  flow*  has  provided  some  insight  into  the  physi¬ 
cal  mechanisms  of  this  instability,  and  rough  information  on  flow  velocity  and  despin  moment.  For 
sufficiently  low  Reynolds  numbers,  more  detailed  results  for  the  velocity  field  have  been  obtained  using 
computational  methods  for  steady  flows6,5.  The  flow  phenomena  at  higher  Reynolds  numbers,  however,  are 
outside  the  scope  of  these  methods,  and  it  is  not  even  clear  whether  the  steady  approach  is  justified. 

Previous  experiments  at  CRDC  and  BRL  were  carried  out  under  full-scale  conditions.  These  studies 
concentrated  largely  on  global  properties  such  as  the  moments  exerted  by  the  internal  fluid  motion.  The 
yet  most  successful  study  of  the  field  properties  is  Miller’s  observation  of  the  void  in  a  partially  filled 
cylinder'.  This  study  shows  ar  axisymmetric  void  at  low  Reynolds  numbers,  a  characteristic  wavy  distor¬ 
tion  of  the  void  at  medium  Reynolds  numbers  and  an  irregular  (probably  unsteady)  liquid-air  interface  at 
high  Reynolds  number.  Computational  studies6  indicate  a  cellular  structure  of  the  flow  at  a  Reynolds 
number  Re  =  45,  where  Re  —  u>a2/u  is  formed  with  the  spin  rate  w,  the  cylinder  radius  a  and  the 
kinematic  viscosity  v.  However,  'here  is  yet  no  link  between  numerical  results  and  void  observations.  An 
attempt  to  trace  buoyant  beads  with  a  movie  camera8  was  very  limited  in  revealing  details  of  the  velo¬ 
city  field.  The  limitations  are  due  to  distortion  of  the  tracer  path  in  the  multi-media  optical  path  involv¬ 
ing  curved  surfaces,  and  to  inevitable  minute  density  differences  in  combination  with  high  accelerations. 
Miller8  used  photochromic  dye  excited  by  a  high-power  pulsed  laser  in  order  to  generate  d  .v.~Grb 
city  profiles.  Lighting  problems  in  recording  the  pictures  by  a  high-speed  movie  camera  forced  a  reduc- 
l.ion  of  the  time  scale,  i.c.  operation  of  the  test  fixture  at  lower  spin  rate,  nutation  rate,  and  kinematic 


"ppendix  A 


in 


viscosity.  Qualitative  pictures  of  the  small  azimuthal  velocity  nave  been  obtained.  The  efforts  to  provide 
more  detailed  data  have  been  discontinued,  however,  due  to  continuing  lighting  problems,  and  the  adverse 
off-design  conditions  at  further  reduced  time  scales. 

In  earlier  work9,  we  have  proposed  a  drastic  reduction  of  length  and  time  scales  for  experimental 
studies,  exploiting  the  principles  of  dynamical  similarity.  Following  these  considerations,  we  have  designed 
and  built  a  low-cost  test  fixture  for  flow  visualization.  In  our  qualitative  approach,  the  length  scale  is 
reduced  to  1/5,  the  time  scale  to  1/10,  thus  reducing  moments  by  more  than  five  orders  of  magnitude  and 
velocities  to  1/50.  In  spite  of  improvising  and  compromising  in  the  interest  of  saving  time  and  money,  we 
have  observed  a  wealth  of  phenomena  from  laminar,  dominantly  unidirectional  flow  through  various 
stages  of  cellular  motions  to  turbulent  motions  with  a  superposed  cellular  structure. 

In  the  following  we  describe  the  principles  underlying  the  design,  the  test  fixture,  the  visualization 
techpiqu*',  find  sorrm  of  our  observations. 


2.  Dimensional  Analysis 

Evaluation  of  the  experimental  attempts  to  visualize  the  fluid  flow  clearly  reveals  the  extreme  full- 
scale  conditions  as  evil.  However,  conclusive  experiments  can  be  conducted  by  exploiting  the  principles  of 
dynamical  similarity  and  appropriate  scaling  laws9'4.  Between  the  three  reference  quantities,  radius  a  , 
spin  rate  u\  and  density  p  for  length,  time,  and  mass,  respectively,  the  density  of  different  fluids  offers  iit- 
tie  variability.  However,  length  scale  and  time  scale  can  be  easily  changed.  For  dynamical  similarity,  the 
lollowing  dimensionless  quantities  must  be  kept  fixed: 


=  e  /'a 

Q 

T  =  0  /u > 

Re  ~  p  uj  a  2//r 


aspect  ratio 
nutation  angle 
frequency 
Reynolds  number 


The  nutation  angle  must  remain  the  same  in  a  scaled  setup.  Radius  a  and  half-length  c  of  the  cylinder 
must  be  scaled  by  the  same  factor  in  order  to  keep  the  aspect  ratio  fixed.  A  second  factor  can  be  applied 
to  both  spin  rate  w  and  nutation  rate  Q  ,  in  order  to  preserve  the  frequency.  Keeping  Re  fixed  requires 
changing  the  kinematic  viscosity  u  =  p  /p  by  the  same  factor  as  tv  a2.  Since  the  desired  tendency  is 
toward  smaller  radii  and  spin  rates,  we  require  less  viscous  fluids  than  those  used  in  the  full-scale  experi¬ 
ments.  Such  fluids  are  easy  to  find. 

It  is  obvious  that  the  main  thrust  of  an  experiment  may  require  specific  optimum  conditions.  Flow 
visualization  requires  low  velocities,  i.c.  low  values  of  iva  .  Measurements  of  moments  require  optimum 
values  of  w2a6.  Minimizing  the  rate  of  change  of  temperature  requires  a  minimum  of  w3a  2.  A  good  setup 
for  flow  visualization,  therefore,  may  produce  moments  in  a  hardly  measurable  range. 

3.  The  Test  Fixture 

The  goal  of  our  efforts  was  to  show  that  a  low-cost  device  (as  $500)  can  be  designed  for  flow  visuali¬ 
zation.  Details  had  to  be  kept  simple.  Accuracy  and  convenience  had  to  compromise.  Various  prelim¬ 
inary  concepts  have  been  condensed  into  the  design  of  a  small  apparatus  that  was  built  and  explored  as  a 
senior  student  project10.  The  result  of  these  efforts  is  sho'  n  in  figure  I.  A  one-inch  inner  diameter 
cylinder  of  aspect  ratio  X  =  4.3  is  used.  The  cylinder  is  cut  from  a  pyrex  glass  tube  with  the  inner  diam¬ 
eter  accurate  within  1/5000  inch,  but  with  varying  wall  thickness  that  affects  the  optical  quality.  1  he 
cylinder  is  filled  with  mixtures  of  water  and  glycerin.  The  mixing  ratio  is  used  to  vary  viscosity.  On  top, 
the  cylinder  is  closed  with  a  screwed-in  plastic  plug.  A  center  hole  allows  access  to  the  interior,  especially 
for  removing  air  bubbles.  The  hole  can  be  closed  using  a  toothpick. 

1  he  cylinder  is  glued  to  a  drive  plug  and  axis  machined  from  a  single  piece  of  aluminum.  The  one¬ 
sided  support  ahows  easy  (optical)  access  to  the  cylinder  and  permits  using  cylinders  of  different  length. 
One-sided  support  is  affordable  due  to  the  moments  being  approximately  five  orders  of  magnitude  smaller 
than  in  the  full-scale  experiments.  The  axis  is  twice  supported  by  ball  bearings.  The  cylinder  and  shaft 
are  driven  via  timing  belts  over  exchangeable  sets  of  pulleys  by  a  <24V  d.c.  motor  with  sufficient  torque 
Hi  the  range  of  500  -  5000  rprn.  Motor  and  cylinder  support  are  mounted  to  an  aluminum  frame  tint  can 
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be  inclined  to  the  vertical  axis  by  approximately  5,  10,  15  and  20*  using  different  support  holes  and 
struts. 

The  horizontal  support  plate  is  machined  to  leave  the  center  position  free  for  access  and  is  screwed 
to  a  commercial  record  player  (Garrard  model  775)  The  plate  can  be  offset  in  order  to  align  the  liquid’s 
center  of  mass  with  the  nutation  axis.  The  record  player  provides  nutation  rates  of  33,  45,  and  78  rpm 
The  hollow  axis  is  utilized  to  provide  power  to  the  spin  motor.  A  nail  with  a  smooth  top  and  a  brush  fixed 
to  the  turntable  proved  sufficient  for  transmitting  a  single  voltage  to  the  motor.  The  remaining  com¬ 
ponents  of  the  experiment  are:  a  Hcathkit-  regulated  power  supply  for  the  spin  motor,  a  strobelight  for 
controlled  pulsed  lighting,  and  suitable  flow  tracers.  The  strobelight  (General  Radio  Strobotac)  with  adju¬ 
stable  frequency  is  used  for  lighting  as  well  as  for  measuring  the  spin  rate  of  the  cylinder. 

4.  Visualization 

As  flow  tracers  we  use  Affiair  100  Silver  Pearl,  kindly  donated  by  EM  Chemicals,  Hawthorne,  NY. 
The  material  consists  of  very  fine  and  shiny  plastic  platelets  commercially  used  for  cosmetic  purposes. 
Although  their  specific  weight  is  different  from  that  of  the  fluid,  the  low  accelerations  in  the  scale  model 
permit  practically  buoyant  behavior  of  the  platelets  over  considerable  time. 

At  the  slow  time  scale  of  the  experiment,  the  fluid  motion  can  be  visually  inspected  while  running 
the  apparatus.  At  high  viscosities,  the  apparatus  can  also  be  suddenly  stopped,  with  the  flow  tracers 
'  frozen”  in  the  resting  fluid.  The  platelets  align  with  surfaces  of  constant  shear.  Therefore,  by  manually 
rotating  the  cylinder  forth  and  back,  the  three-dimensional  structure  of  the  field  can  be  inspected.  This 
crude  observation  is  very  helpful  in  developing  the  visualization  technique.  A  detailed  account  of  the 
technique  (appropriate  particle  density,  pitfalls  such  as  the  history  of  particle  distribution  and  alignment) 
has  been  given  elsewhere10. 

Visualization  of  the  frozen  pattern  can  be  essentially  improved  by  using  a  light  sheet  passing 
through  the  spin  axis.  Sheet  lighting  enhances  the  clarity  of  the  flow  pattern  by  showing  only  the 
reflecting  particles  in  a  cut  through  the  fluid.  It  reduces  the  undesirable  reflections  from  the  cylindrical 
surfaces  and  also  enables  photographic  recording  of  the  flow  structure  while  the  apparatus  is  in  operation. 
A  continuous  light  sheet  is  produced  by  a  Spectra  Physics  model  120  (15  mWl  helium-neon  laser  and  a 
cylinder  lens.  In  order  to  avoid  the  need  for  accurately  firing  the  camera  (3a  mm  Pentax  with  50  mm 
lens)  at  a  certain  time,  a  cylindrical  card  board  screen  with  a  vertical  slot  and  a  90*  offset  opening  is  fixed 
to  the  circumference  of  the  turntable.  The  shutter  is  manually  opened  and  closed  after  the  laser  sheet  of 
light  Hashed  3  to  5  times  through  the  slot. 

5.  Results 

Some  photographs  taken  with  the  apparatus  in  motion  are  shown  in  Figures  2-7.  The  figures  show 
the  flow  pattern  in  the  plane  spanned  by  spin  axis  and  nutation  axis  for  0  =  21.3* ,  ft  =-78  rpm  and 
different  Reynolds  numbers.  Figure  2  shows  that  at  Reynolds  numbers  as  low  a s  Re  =  20  a  cellular  pat¬ 
tern  develops  with  a  pronounced  symmetry  about  the  axis  as  well  as  the  midplane  of  the  cylinder  At  the 
present  time  it  is  unclear  whether  this  pattern  reflects  the  instantaneous  velocity  field  Symmetry  argu¬ 
ments  support  viewing  this  pattern  as  origin*/ing  from  a  nonlinear  streaming  term  As  Re  increases  to 
Rt  =  40  (Figure  3),  the  pattern  and  its  symmetry  become  more  pronounced.  At  Re  =  50  (Figure  4), 
additional  cells  develop  near  the  cylinder’s  midplane.  Simultaneously,  the  symmetry  with  respect  to  the 
cylinder  axis  is  broken.  A  characteristic  wavy  distortion  of  the  pattern  near  the  axis  develops  that  is 
more  clearly  shown  in  Figure  5  at  Re  =  105.  While  the  cells  disappeared,  virtually  axisymmetric  hub¬ 
bies  occur  at-  the  end  plates.  At  Re  =  140  (Figure  6)  these  bubbles  still  persist.  The  bright,  wavy  line 
near  the  axis  has  broken  into  segments  that  are  very  much  aligned  like  the  void  in  Miller’s  observations. 
This  pattern  occurs  only  in  the  plane  of  spin  axis  and  nutation  axis  and  is  therefore  considered  to 
represent  the  instantaneous  velocity  field.  From  the  wealth  of  increasingly  complex  phenomena,  Figure  7 
finally  shows  a  visualization  at  high  Re  =  8000  The  random  distribution  of  the  particles  in  the  interior 
most  likely  indicates  turbulent  flow.  Nevertheless,  the  faint  line  near  the  axis  resembles  the  characteristic 
centerline  distortion  of  Figure  C,  indicating  a  superposed  large  scale  structure.  The  presence  of  such  a 
large  scale  motion  is  also  supported  by  the  regular  bands  of  particles  deposited  at  the  cylinder  wall 
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Specification  of  accurate  Reynolds  numbers  suffers  from  some  uncertainty  in  monitoring  and  measur¬ 
ing  the  wide  range  of  viscosities  for  the  hygroscopic  water-glycerin  mixtures  exposed  to  uncontrolled  ther¬ 
mal  conditions  To  within  this  uncertainty,  however,  the  figures  clearly  reveal  the  cellular  structure  of  the 
flow  and  the  changes  of  the  structure  as  the  Reynolds  number  increases.  Perhaps  the  most  striking  result 
of  this  visual  otudy  of  the  flow  st.uct-ure  is  the  manifold  of  pattern  at  higher  Reynolds  numbers.  A  sys¬ 
tematic  analysis  of  these  patterns  has  not  been  conducted.  Although  we  found  numerous  opportunities  for 
improvements,  the  feasibility  of  flow  visualization  with  relatively  simple  means  by  proper  scaling  has  been 
clearly  demonstrated. 
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Figure  1 


The  miniature  test  fixture.  The 
inner  radius  oi'  the  cylinder  is 
a  —  1.27  cm,  the  aspect  ratio  4.3. 


Figure  3 

0  —  21. .r 

fi  =  78  rpm 
v  =  240  cSt 
w  =  600  rpm 
Re  =  40 
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Figure  7 

0  —  21.3  * 
ft  =  78  rpin 
u  ■=  0.95  cSt 
w  ==  450  rpm 
Rc  =  8000 
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ABSTRACT 

Spin-stabilized  projectiles  with  liquid  payloads  can  experience  a  severe  flight  instability 
characterized  by  a  rapid  yaw-angle  growth  and  a  simultaneous  loss  in  spin  rate.  Labora¬ 
tory  experiments  and  field  tests  have  shown  that  this  instability  originates  from  the 
internal  fluid  motion  in  the  range  of  small  Reynolds  numbers.  In  earlier  work,  we 
developed  a  simple  model  of  this  flow  based  on  linearized  equations  for  the  deviation 
from  solid-body  rotation  in  an  infinite  cylinder.  Here,  we  perform  a  perturbation 
analysis  in  order  to  estimate  the  effect  of  nonlinear  terms.  Beyond  a  small  correction  of 
the  social  velocity  component,  we  obtain  radial  and  azimuthal  components  of  the  velocity 
field  in  agreement  with  computational  results  for  the  core  region  of  a  finite-length 
cylinder.  The  analytical  results  are  exploited  in  the  design  of  a  spectral  Navier-Stokes 
solver  for  the  steady  motion  in  a  finite  cylinder.  A  first  raw  version  of  this  spectral  code 
provides  flow  field  and  pressure  distribution  in  a  small  fraction  of  the  computer  time 
required  by  existing  codes.  We  report  some  results  and  discuss  possible  refinements  of 
this  code. 

1.  Introduction 

It  is  well-known  that  spin-stabilized  shells  carrying  liquid  payloads  can  suffer  a 
dynamical  instability  which  results  in  an  increased  coning  (or  yaw)  angle  and  a  simul¬ 
taneous  loss  in  spin  rate  Laboratory  experiments,  computational  results,  and  field  tests 
indicate  that  these  phenomena  arise  from  the  coning-induced  fluid  motion  in  a  limited 
range  of  small  Reynolds  numbers.  Although  in  special  cases  this  instability  has  been 
removed  by  trial  and  error,  future  design  of  reliable  projectiles  would  profit  from  the 
opportunity  to  estimate  the  liquid  moments,  and  to  include  these  moments  in  flight 
simulators.  The  empirical  data  base  [l,  2)  is  sparse,  however,  and  computational 
methods  in  use  (3,  4,  5j  are  r«.  .  sr  demanding. 

Our  theoretical  analysis  of  this  problem  serves  on  one  hand  to  gain  insight  into  the 
anatomy  of  the  flow  phenomena  and  to  support  the  ongoing  experiments.  On  the  other 
hand,  it  promotes  our  efforts  to  develop  a  more  efficient  code  for  the  numerical  simula¬ 
tion  of  the  flow  in  a  fini'  container.  While  the  analytical  work  aims  at  the  velocity  field 
in  the  core  region  of  a  sufficiently  long  cylinder  and  on  the  viscous  components  of  the 
moments,  in  particular  the  viscous  despin  (negative  roll)  moment,  the  computational 
vork  also  captures  the  flow  near  the  end  walls  and  the  pressure  contributions  to  yaw 
and  pitch  moments. 

Our  previous  work  [6]  shows  that  the  deviation  from  solid  body  rotation  is 
governed  by  a  small  parameter  (  —  f2.sin0/u/  involving  the  nutation  rate  fi  ,  the  nuta¬ 
tion  angle  0 ,  and  the  spin  rate  w.  The  solution  of  the  linearized  equations  consists  of 
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only  an  axial  component  of  order  0(e).  This  axial  flow  is  the  dominating  feature  of  the 
fluid  motion  and  produces  a  negative  roll  moment  of  order  0(e2)  owing  to  Coriolis 
forces.  Although  these  results  are  in  reasonable  agreement  with  experimental  and  com¬ 
putational  data,  one  may  anticipate  modifications  of  velocity  field  and  roll  moment  if 
nonlinear  terms  are  taken  into  account.  Estimates  of  these  nonlinear  effects  are  desired 
in  order  to  support  previous  results  and  to  verify  our  conclusion  that  the  three- 
dimensional  flow  field  in  a  finite-length  cylinder  is  essentially  given  by  the  solution  of 
linearized  momentum  equations.  In  the  following,  we  perform  a  straightforward  pertur¬ 
bation  expansion  for  the  nonlinear  problem.  We  develop  and  solve  the  equations  for  the 
flow  in  an  infinitely  long  cylinder  up  to  order  0(t3).  A  closed-form  solution  is  given  for 
the  radial  and  azimuthal  velocity  components  at  second  order.  The  third-order  equa¬ 
tions  arc  solved  numerically. 

The  perturbation  solution  also  provides  estimates  for  the  number  of  expansion 
functions  required  for  accurate  spectral  representation  of  the  radial  (r)  and  azimuthal 
(0)  structure  of  the  solution.  A  spectral  code  appears  as  an  attractive  alternative  to  the 
existing  Navier-Stokcs  solvers.  The  finite-difference  code  developed  at  Sandia  Labora¬ 
tories  [3,  1]  exploits  Chorin’s  method  of  artificial  compressibility.  The  steady  solution  at 
11  X  2d  X21  grid  points  in  r  ,  <t> ,  z -direction  is  obtained  by  integrating  over  typically 
l O4  time  steps,  a  task  that  requires  68  minutes  of  CPU  time  on  an  IBM  3090.  The  result 
consists  of  22,000  plus  values  for  the  velocity  components  vr ,  v+,  vt  and  the  pressure  p 
that  can  be  utilized  for  a  calculation  of  the  moments.  Strikwerda  &  Nagel  [5]  describe  a 
code  using  finite  differences  in  radial  and  axial  direction  and  pseudospectral  differencing 
in  the  azimuthal  direction.  Nonuniform  grids  are  introduced  for  increased  resolution 
near  the  walls.  The  difference  equations  are  solved  by  an  iterative  method  based  on  suc¬ 
cessive  over-relaxation.  The  computer  time  required  is  comparable  to  that  of  the  Sandia 
code  (Nusca,  BRL,  personal  communication).  Although  the  relative  merits  of  the  two 
codes,  especially  with  respect  to  the  captured  range  of  Reynolds  numbers  are  yet  in  the 
dark,  it  seems  well  possible  to  beat  both  of  these  codes  in  two  respects:  computer  time 
and  adaptability  to  the  unsteady  problem. 

For  a  feasibility  study,  we  have  pursued  a  simple  concept  that  is  open  to  numerous 
refinements.  We  use  Ohebyshov-Fourier-Chebyshev  expansions  in  r,d>,z,  respectively, 
and  convert  the  linearized  equations  into  a  linear  algebraic  system  for  the  expansion 
coefficients  The  solution  of  this  system  (or  any  other  solution  at  neighboring  parame¬ 
ters)  is  used  as  initial  approximation  for  iterative  improvement  by  the  modified  Newton 
method.  The  experience  with  this  co-.e  is  encouraging  with  respect  to  accuracy, 
efficiency,  and  robustness. 

2.  Governing  Equations 

We  consider  the  motion  of  a  fluid  of  density  p  and  viscosity  p  in  a  cylinder  of 
radius  a  and  length  2c  that  rotates  with  the  spin  rate  u  about  its  axis  of  symmetry,  the 
c-axis.  We  consider  the  motion  with  respect  to  the  nutating  coordinate  system  z  ,  y,  z. 
This  system  is  obtained  from  the  inertial  system  X .  Y ,  Z  by  a  rotation  with  the  nuta¬ 
tion  angle  0  about  the  axis  Y  y  Therefore,  x  is  in  the  Z ,  z -plane,  and  this  plane 
rotates  about,  the  Z-axis  with  the  nutation  rate  H  .  The  two  axes  of  rotation  intersect 
in  the  center  of  mass  of  the  cylinder  We  consider  >  0,  fi  ,  and  0  <  0  <  fr /2  as  con- 
>t  :mt 


a 

d-> 
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The  fluid  motion  is  governed  hv  the  Navier-Stokes  equations  written  in  the  nutat¬ 
ing  coordinate  system: 

+  x  Vn  +  n  X  (n  X  r)j  =  -  ypn  i-mvX  .  (la) 

v  v„  ----  o  (lb) 

V„  is  the  velocity  measured  in  the  nutating  frame,  Pn  the  pressure,  and  r  the  position 
vector.  Equations  (1)  are  subject  to  the  no-slip  and  no- penetration  conditions  at  the 
cylinder  walls. 

It  is  convenient  |6)  to  split  the  velocity  and  pressure  fields  according  to 

Vn  -  V,  -r  V,  .  l\  Pt  h  Pd  ,  (2) 

where  V,,  P,  describe  the  state  of  pure  solid-body  rotation,  whereas  V^,  Pd  represent 
the  deviation  from  solid- body  rotation.  The  deviation  and  the  reduced  pressure  Pd 
are  ultimately  responsible  for  the  observed  flight  instability 

The  equations  for  V^,  Pd  are  written  in  terms  of  nondimensional  quantities  ,  pd 
using  o  ,  w,  and  p  for  scaling  length,  time,  and  mass,  respectively.  The  solution  then 
depends  on  four  nondimensional  parameters:  aspect  ratio  X  ---  c  !a,  nutation  angle  6, 
frequency  r  —  U  /co,  and  Reynolds  number  R  —  puia'/p.  The  aspect  ratio  enters  the 
solution  only  through  the  boundary  conditions  at  the  end  walls  of  the  cylinder  The 
boundary  conditions  on  are  homogeneous. 

In  cylindrical  coordinates  r  ,  6 .  z  ,  the  equations  for  the  nondimensional  deviation 
velocity  \d  =  (i>r  ,  v, )  and  pressure  pd  take  the  form 

1  d  ,  w  1  dv*  x  dvt 

— 3— (rvr )  -b - rr-  +  -7-  =  0  ,  (3a) 

r  dr  r  a<p  az 

v  1 

D'vr  -  ——  2(1  f  rt) v p  f  2rd>vl  (3b) 

dPi  ,  1  iM„  ‘V  2  dv* 

-  -7—  -r  -=\D  vr  -  -T-——  j  , 

or  n  ri  r*  o <*> 

D'v.  ■+•  — - — —  -i  2(1  -I  r  )vr  -  2rrv:  (3c) 

r 


r  86  R  ' 


!±.  ■  jl  ,)v’  1 

r2  ^  r2'  **  1  ' 


D'v,  -1-2 T rV t  -  2 TjV,  ■- 


2  rr,  +  ±D"  vt 


w  ne  re 


'() 

(J 

8 

-J  t 1  - 

v*  (J 

8 

\  _ 

(Jl 

8  6 

1  lr  r 

(Jr 

r  86 

’’  iJz 

U 2 

1  (1 
r  (Jr  ^ 

1  d‘  3 

8l 

(It  2 

2  '1  1  2 

r  i  86 

fJz 2 

Appendix  ll 


20 


(4) 


rr  =  -ecos <f>  ,  Tf  ~  £sin<^  ,  t,  =  tcosO  ,  £  =  rsin#  . 

The  primary  effect  of  nutation  is  contained  in  the  ^-periodic  force  term  -2rrr  = 
2t  r  cos<t>  in  the  2-momentum  equation  (3d).  For  e  =  0,  equations  (3)  support  the 
trivial  solution  vd  =  0 ,  pd  =  0.  The  system  also  supports  the  following  symmetries: 


tV (r  ,  z  )  —  vr{r  ,4>,  z) 

(5a) 

v  ^(r  ,  <t>+ir ,  -  z)  =  Vf{r  ,<t>,  z) 

(5b) 

vi  (r  ,  <t>+n,  -  z)  =  -  vt(r  ,  <f>,  z) 

(5c) 

Pi{r  ,  4>+tt,  -  z)  =  pj(r  ,<t>,  z) 

(5d) 

3.  Perturbation  analysis  for  an  infinite  cylinder 

The  steady  flow  in  a  relatively  long  cylinder  (aspect  ratio  X  >  4)  at  low  Reynolds 
number  is  expected  to  have  a  rather  simple  structure  and  to  provide  a  roll  moment  pro¬ 
portional  to  Re  .  In  fact,  the  flow  is  expected  to  exhibit  little  axial  variation  over  much 
of  the  cylinder  length  Previous  work  (Gj  has  therefore  relaxed  the  boundary  conditions 
at  the  end  walls.  In  this  way,  one  seeks  the  steady  flow  in  a  finite  segment  of  an 
infinitely  long  cylinder. 

In  the  physical  situations  of  interest,  £  —  (H  /u>)  sin0  is  a  small  parameter, 
(  <  0  06  Consequently,  it  seems  reasonable  to  pursue  a  straightforward  perturbation 
expansion  in  c  This  provides  in  the  form 

v,  =  £  ^{n\r  ,4>)  (6) 

n  —  1 

and  similar  expressions  for  pd  . 

The  development  of  general  expressions  for  the  expansion  coefficients  from 
equations  (3)  indicates  an  alternating  pattern:  Odd-order  terms  contain  odd  multiples  of 
<!>  and  contribute  only  to  the  axial  velocity  vt ,  while  even-order  terms  contain  even  mul¬ 
tiples  of  the  azimuthal  coordinate  <f>  and  contribute  only  to  the  radial  velocity  vr  and 
azimuthal  velocity  v^.  Therefore, 


(0,  0,  t//n^),  n  odd  , 

Kln),  V#  {n  \  0),  71  even  , 


and  the  components  of  v'n  *  take  the  form 


«V( "  >  -  "s  (v(r)e’2^  4-  dnm(r  )  e-,2m*)  , 

m  --- 1 


n  /2 


tV"1-  vno{r)  +  £  (vnm  (r  )  e>2m*  4-  vnm  (r  )  e  '  Zm  *)  , 

m  1 

HI  1 


(7) 


(8a) 

(8b) 

(8c) 


where  the  tilde  denotes  the  complex  conjugate.  The  aperiodic  term  in  is  suppressed 
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by  the  continuity  equation.  The  r -dependent  coefficient  functions  in  eqs.  (8)  are 
required  to  satisfy  homogeneous  boundary  conditions  at  r  —  1  and  to  be  finite  at  the 
axis  r  =  0  for  a  physically  meaningful  solution. 


At  the  lowest  order  0{t),  the  * -independent  force  term  in  eq.  (3d)  can  be  balanced 
only  by  an  axial  component  of  the  deviation  velocity.  This  component  is  the  dominating 
feature  of  the  flow  in  a  long  cylinder.  The  axial  velocity  at  order  0(e)  can  be  found  in 
analytical  form, 

i(-w~ r)’  (10) 

where  is  the  modified  Bessel  function,  and  o  =  (1  +  i  )(R  j  2)1,/2.  This  solution  is 
valid  for  arbitrary  Reynolds  number  but  may  be  unstable  as  R  exceeds  some  critical 
value.  The  properties  of  .he  resulting  flow  field  are  discussed  by  Herbert  [6]. 

At  higher  order,  it  is  convenient  to  eliminate  the  pressure  for  the  periodic  com¬ 
ponents  by  using  the  vorticity  form  of  eqs.  (3).  At  order  0(e2),  comparison  of  the  equa¬ 
tion  for  t/30  with  the  imaginary  part  of  the  equation  for  wn  immediately  shows  that  the 
aperiodic  component  of  the  azimuthal  velocity  is 


VvJir)  —  -2  Im(ien(r  )).  (11) 

This  relation  can  be  exploited  to  show  that  the  despin  moment  of  order  0(<2)  due  to 
shear  forces  on  the  cylinder  wall  is  identical  with  our  former  result. 

The  ^-periodic  components  are  governed  by  a  coupled  set  of  inhomogeneous 
differential  equations  with  variable  coefficients.  Essential  simplification  at  the  expense  of 
increasing  the  order  of  differentiation  results  from  eliminating  t>2l  by  use  of  the  con¬ 
tinuity  equation.  With  some  effort,  the  radial  velocity  component  of  0(e2)  can  be  found 
in  closed  form, 


u2i(a)  =  -[c^s)  +  c^a)]  +  c3a  +  ~-  4-  ^ 

a  a*  a 


Ji(0/V2) 


(12) 


where  a  —  0r  ,  0  =  (i  -  l)/?1/2,  and  J  2,  andV2  are  Bessel  functions.  The 
coefficients  c  i,  c2,  c3,  and  c4  can  be  determined  numerically. 

The  velocity  components  at  order  0(e3)  are  of  interest  primarily  since  w3l  provides 
the  first  nonlinear  correction  to  the  despin  moment.  In  view  of  the  effort  involved  in 
deriving  the  closed  form  solution  for  w21  and  the  ultimate  need  to  determine  the 
coefficients  in  eq.  ^12)  numerically,  we  decided  to  solve  the  differential  equations  for  the 
third-order  components  by  means  of  a  spectral  collocation  method. 


4.  Results  of  tiie  Perturbation  Aoialysis 

Detailed  equations,  results,  and  graphs  of  the  various  functions  at  relevant  Rey¬ 
nolds  numbers  will  be  published  elsewhere  [7].  Here  we  give  only  a  summary  of  the  main 
results.  The  motion  is  governed  by  the  axial  component  u/ u  at  order  0(e)  Of  the 
higher  order  terms,  only  the  aperiodic  term  is  substantial.  In  the  cylinders  center 
section,  these  terms  are  in  gooo  agreement  with  results  obtained  from  the  Sandia  code, 
and  in  excellent  agreement  with  our  own  computations.  All  the  other  terms  are  not  only 
of  order  0(1)  but  in  fact  less  than  unity,  assuring  rapid  convergence  of  the  perturbation 
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series.  The  contribution  of  tu31  to  the  despin  moment  is  negligible.  The  ^-periodic 
terms  oscillate  about  zero  as  r  varies  between  0  <  r  <  1.  Accurate  representation  of 
single  high-order  terms  by  radial  Chebyshev  series  may  require  numerous  expansion 
functions  For  the  total  velocity  field,  however,  the  error  in  representing  these  terms  is 
of  little  importance.  At  Reynolds  numbers  in  the  range  of  maximum  despin  moment, 
reasonably  accurate  approximations  can  be  obtained  with  as  few  as  five  polynomials  in 
radial  direction.  In  the  azimuthal  direction,  the  solution  is  governed  by  terms  periodic  in 
<?,  and  by  the  aperiodic  term  v20.  Fourier  series  with  three  or  five  modes,  therefore,  pro¬ 
vide  approximations  of  sufficient  accuracy  for  practical  purpose. 

5.  Spectral  Approximations  for  a  Finite-Length  Cylinder 

The  results  of  the  perturbation  analysis  suggest  that  a  good  approximation  to  the 
flow  in  a  finite  cylinder  can  be  obtained  by  solving  the  linearized  version  of  equations  (3). 
Linearization  can  be  performed  in  different  ways.  The  first  is  a  linearization  in  e,  as  in 
the  perturbation  analysis.  The  resulting  equations  support  strong  symmetries.  Beyond 
equations  (5),  the  solution  satisfies 

vd{r  ,  <j>+n ,  z)  =  -vd(r  ,  <j>,  z)  ,  (13a) 

Pi{r  ,  4>+n,  z)  =  -Pi(r,<f>.  z)  .  (13b) 

These  relations  provide  a  useful  check  on  the  results  of  the  spectral  code.  A  second 
linear  system  can  be  obtained  by  linearization  in  the  components  of  vd  .  This  lineariza¬ 
tion  retains  coupling  terms  such  as  2t  ^v,  in  eq.  (3b)  which  destroy  the  symmetries  (13). 
The  second  system  can  be  considered  a  special  case  of  a  third  linearization  about  some 
known  solution  v]°\  pj°l  The  third  procedure  is  very  efficient  if  the  solution  is  sought 
for  a  densely  spaced  sequence  of  parameter  combinations  as  in  flight  simulations.  The 
second  system  is  equivalent  with  the  third  one  for  vj°)  ==  pd ^  =  0. 

The  algebraic  form  of  the  equations  is  obtained  by  use  of  spectral  collocation.  The 
velocity  components  are  expressed  in  the  form 

K  L  M 

t V=£  £  £  uklmRZ{r)FM)Zm{z/\)  (14) 

k  |  l  =»1  m  =  1 

with  similar  expressions  for  v0,  vt,  and  pd .  The  azimuthal  functions  are 
/'i  --  cos  '[(/  -  \)4> /2\  for  odd  l ,  Ft  —  sin  (/  <^» /2j  for  even  /,  where  1  =  1,  2,  •  •  •  L  , 
and  L  is  odd  The  azimuthal  collocation  points  are  equidistant,  4>t  =  2zr(/  -  1  )/L  .  If 
no  use  of  the  symmetries  (5)  is  made,  the  axial  expansion  functions  are  the  Chebyshev 
polynomials  Zm  —  Tm  _  j(z  /X),  m  =  1,  2,  •  •  ■  M.  The  collocation  points  are 
zm  X  —  cos  j(m  -  1)t  j{M  -  1)|.  In  radial  direction,  even  or  odd  Chebyshev  polynomi¬ 
als  are  used,  depending  on  the  quantity  under  consideration  and  the  periodicity  in  <f> . 
The  proper  choice  is  dictated  by  the  requirement  of  a  unique  value  of  all  quantities  on 
the  axis  r  —  0  For  example,  the  axial  velocity  component  must  assume  a  unique  value 
independent  of  <*>  as  r  — ►  0.  Therefore,  even  polynomials  are  to  be  used  if  /  =  1  while 
odd  polynomials  are  to  be  used  it  l  >1  The  radial  collocation  points  are  r*  = 
cos  \(k  -  1 ) rr /( 2 K  -  1)|,  k  —  1,2,  K  Consequently,  0  <  rk  <1,  and  no 
difficulty  can  arise  from  points  on  the  axis.  The  collocation  points  in  radial  and  axial 
direction  are  concentrated  near  the  boundary  such  that  high  resolution  in  this  region  is 
obtained  without  additional  coordinate  transformations. 
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Our  implementation  of  the  spectral  method  uses  precalculated  and  stored  matrices 
containing  the  values  of  the  expansion  functions  and  their  derivatives  at  the  collocation 
points.  It  is  a  straightforward  matter  to  convert  the  linear  system  of  part  ml  differential 
equations  derived  from  eqs.  (3)  into  an  algebraic  system  of  dimension  A  —4  K  L  M 
for  the  coefficients  uktm  ,  vklm,  wklm  ,  and  pki,n  for  vr  ,  v^,  vz  ,  and  pd  ,  respectively.  It  is 
not  straightforward,  however,  to  implement  the  homogeneous  boundary  conditions  for 
the  velocities  at  the  cylinder  wall  and  the  condition  on  the  pressure  that  is  only  deter¬ 
mined  to  within  an  additive  constant.  In  principle,  the  boundary  conditions  are  imple¬ 
mented  by  replacing  three  of  the  four  differentia!  equations  in  the  boundary  points.  The 
question  then  is  which  equation  should  be  retained  and  where  the  condition  on  the  pres¬ 
sure,  e  g  pd  —  0,  should  be  applied.  Trial-and-etror  leads  to  numerous  cases  with  ill- 
determined  matrices  or  zero  determinant.  In  other  cases,  a  correct  solution  for  the  velo¬ 
city  field  is  obtained,  but  the  pressure  contains  a  non-physical  spurious  term  With  the 
velocity  field  given,  we  attempted  to  calculate  the  pressure  by  solving  a  Poisson  equation 
with  von  Neumann  boundary  conditions,  but  we  encountered  the  same  difficulties. 
Problems  with  calculating  the  pressure  in  closed  domains  with  spectral  methods  are 
well-known,  e  g.  j8j.  However,  tiie  reports  of  negative  results  are  rather  unspecific,  and 
neither  the  origin  nor  methods  for  removal  of  this  spurious  term  seem  to  be  known. 

We  have  therefore  performed  a  detailed  analysis  of  the  flow  in  a  square  driven  by 
an  internal  force  field  This  simpler  two-dimensional  problem  exhibits  all  characteristics 
-  including  the  spurious  pressure  term  -  of  the  original  problem  Detailed  results  of  this 
study  will  be  reported  elsewhere  [9j.  The  study  reveals  that,  the  spurious  term  is  associ¬ 
ated  with  the  corners  of  the  domain.  The  term  vanishes  in  all  collocation  points  except 
the  corners,  where  it  may  assume  arbitrary  values  The  term  can  be  suppressed  by 
retaining  in  the  corners  one  of  the  momentum  equations  that  contain  the  derivative  of 
the  pressure  in  the  direction  of  the  boundary.  In  the  cylinder  problem,  the  2- 
momentum  must  be  retained  in  order  to  suppress  even  as  well  as  odd  spurious  terms. 
The  condition  on  the  pressure  can  be  applied  anywhere  except  in  the  corner  points. 

We  solve  the  linear  algebraic  system  for  the  expansion  coefficients  with  a  special 
subroutine  based  on  Gauss  elimination  with  partial  pivoting  The  subroutine  stores  all 
data  required  to  solve  the  same  system  with  a  new  right-hand  side  without  repeating  the 
costly  (0  (yV3)  operations)  reduction  of  the  matrix  to  upper  triangular  form.  Once  the 
solution  is  obtained,  a  new  right-hand  side  is  formed  taking  the  nonlinear  terms  into 
account  and  the  system  is  solved  again.  This  procedure  is  iteratively  repeated  until 
sufficient  accuracy  is  obtained.  The  procedure  is  equivalent  to  the  modified  Newton 
iteration  (without  updating  the  Jacobian  in  every  step)  and  converges  rapidly  since  the 
nonlinear  corrections  to  the  velocity  are  small  while  the  pressure  appears  linear  in  equa¬ 
tions  (3). 

0.  Results  oi  the  Spectral  Code 

In  the  following,  we  present  some  preliminary  results  of  a  test  run  for  H  =  14  95, 

0  20c  ,  r  --  O.lGfiT.  and  X  --  4.3(>8  which  results  in  t  —  0.057.  The  results  are  for 
K  =  4,  L  —  ,\f  —  5.  and  consequently  \  =-  400  Detailed  convergence  tests  will  be 
performed  with  later  versions  of  the  spectral  code.  Figure  1  shows  the  axial  and  radial 
velocity  in  the  x,  2-plane.  Only  the  upper  half,  z  >  0,  of  the  cylinder  is  shown;  the 
lower  half  is  governed  by  the  symmetries  (5).  The  velocity  distribution  at  2  =  0  agrees 
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well  with  the  results  of  the  perturbation  analysis  and  computations  with  the  Sandia 
code.  Near  the  walls,  the  solution  seems  to  be  more  realistic  and  more  accurate  than  the 
Sandia  results  The  figure  also  verifies  the  existence  of  a  predominantly  axial  flow  over 
most  of  the  cylinder  length,  except  within  a  region  of  the  order  of  the  radius  near  the 
end  wall.  Linear  and  nonlinear  velocity  distributions  are  hardly  distinguishable.  Clearly 
visible  is  the  turning  of  the  flow  near  the  end  wall.  The  radial  and  azimuthal  velocities 
at  2  -•=  0  9\  are  shown  in  figure  2.  The  right  tick  mark  indicates  the  x -direction, 

0  —  0  At  Re  =  14.95,  the  maximum  of  the  axial  velocity  occurs  at  45 

Pressure  distributions  in  the  x  ,  z -plane  are  given  in  figures  3  and  4  with  the  heavy 
lines  indicating  positive  values.  The  pressure  in  figure  3  is  obtained  simultaneously  with 
v i  from  equations  linearized  in  c,  and  clearly  shows  the  symmetry  (13b).  Figure  4  gives 
the  result  from  the  nonlinear  equations.  It  is  interesting  to  note  that  a  very  similar  pres¬ 
sure  field  can  be  obtained  by  solving  the  Poisson  equation  for  the  pressure  with  the 
linear  velocity  field.  The  inhomogeneous  term  in  the  Poisson  equation  is  inherently  non¬ 
linear  in  the  velocities.  Figure  5  gives  the  pressure  distribution  across  the  cylinder  near 
the  end  wall  at  ^  =  0  9\  Remarkable  is  the  formation  of  a  high-pressure  region  in  the 
corner  near  o  ---  0,  which  produces  a  large  moment  about  the  y-axis.  Looking  at  a  series 
of  plots  like  figures  4  and  5,  one  may  wonder  whether  the  details  of  the  pressure  varia¬ 
tion  near  the  cylinder  wall  can  be  resolved  with  a  finite  difference  approximation  with  a 
step  size  of  Ar  =0.1. 

The  azimuthal  mean  velocity  at  c  —  0  is  shown  in  figure  6.  The  shear  exerted  by 
this  component  on  the  cylinder  wall  opposes  the  spinning  motion  and  is  the  ultimate 
cause  of  the  despin  moment.  The  axial  and  radial  mean  velocity  field  is  given  in  figure 
7  This  streaming  term  exhibits  a  toroidal  motion  stretched  over  each  half  of  the 
cylinder.  It  is  this  mean  velocity  that  causes  the  symmetric  pattern  in  flow  visualiza¬ 
tions  [101.  Figure  8  shows  the  observed  pattern  of  the  flow  at  R  ^  30  which  is  typical 
lor  the  range  of  low  Reynolds  numbers. 

7.  Discussion 

The  experience  with  the  first  version  of  the  spectral  code  shows  that  high  perfor¬ 
mance  can  be  achieved.  The  reported  run  with  N  —  400  requires  1.3  minutes  CPU 
time  on  an  IBM  3090,  48  minutes  on  an  Apollo  DN300  desktop  computer.  The  solution 
is  obtained  in  semi-analy lical  form  with  only  ;V  =  400  numerical  coefficients.  This  low 
data  volume  is  especially  attractive  for  communication  with  remote  supercomputers. 
The  code  is  very  well  suited  for  vectorization,  since  practically  all  CPU  time  is  spent  on 
constructing  and  solving  an  algebraic  system.  However,  the  code  demands  larger 
memory  than  other  codes  [3,  5|.  Since  64-bit  arithmetic  is  highly  recommended  for  spec¬ 
tral  methods  in  general,  and  the  algebraic  system  requires  ,V(,V  f  l)  words  of  storage, 
the  above  test  requires  13  Mbyte  of  memory  Nowadays,  the  memory  requirement 
appears  acceptable  even  if  higher  resolution  is  desired 

Finally,  there  are  various  ways  to  improve  the  performance  and  lessen  the 
demands.  The  first  step  is  to  exploit  symmetry  which  reduces  ,V  by  a  factor  of  1/2, 
storage  by  11.  and  time  by  =s  18  Second,  the  solution  process  can  be  split  into  two 
levels,  the  first,  of  which  calculates  only  the  velocity  components  while  the  pressure  is 
obtained  a  posteriori  by  solving  the  Poisson  equa* ion.  After  these  changes,  the  above 
lost  run  w  e|  ,nlniiT  less  i.han  I  miniile  on  an  \1( '68020  68881  based  desktop  computer 


Appendix  r 


25 


Alternatively,  runs  with  higher  resolution  can  be  executed  within  a  short  time  on  super¬ 
computers.  One  may  also  consider  reducing  the  storage  requirement  by  line  iteration. 
However,  the  ability  to  obtain  a  reasonably  acur&te  solution  by  direct  solution  of  the 
(large)  algebraic  system  bears  valuable  potential  to  answer  the  question  whether  the 
steady  solution  is  stable,  and  allows  for  analysis  of  unsteady  motions.  Tv.e  design  of  a 
reliable  code  for  the  unsteady  problem  can  take  profit  from  the  kowledge  of  the  eigen¬ 
value  spectrum  for  small  unsteady  disturbances  of  the  steady  flow. 
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Figure  Vector  plot  of  the  velocity  field  Figure  2.  Vector  plot  of  the  velocity  field 

the  i ,  c  -plane  for  c  >  0.  across  the  cylinder  at  2  /X  -  0.9. 
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Figure  4.  Contour  plot,  of  t.lic  nonlinear 
pressure  Held  in  the  1,2-plane  for 
2  >  0 
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Figure  5.  Contour  plot  of  the  pressure 
field  across  the  cylinder  at  z  /X  =  0.9 


Figure  6.  Vector  plot  of  the  mean  velo¬ 
city  field  across  the  cylinder  at  z  /\  —  0. 


Figure  7. 
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Vector  plot  of  the  mean  veto-  Figure  8.  Pattern  of  the  fluid  motion  at 

in  the  z  ,  z -plane  for  z  >  0.  low  Reynolds  numbers  (/?  30)  in  the 

x  ,  z  -plane. 
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ABSTRACT  • 
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■4 

Artillery  shells  with  liquid  payloads  may  experience  a  severe  (light  instability  owing  to  the  viscous  fluid 
motion  in  the  cylindrical  payload  container.  Analytical  studies  of  this  fluid  motion  suggest  that  the  j 
velocity  field  in  the  aeroballistic  coordinate  system  can  be  obtained  from  a  linearized  system  of  cqua-  i 
tions,  with  a  small  correction  for  nonlinear  terms.  Moreover,  the  smooth  solution  of  the  Navier-Stokes  j 
equations  at  the  relevant  low  Reynolds  numbers  can  be  accurately  approximated  by  relatively  few  | 
terms  of  spectral  expansions.  We  describe  a  spectral  collocation  method  for  calculating  velocity  and  I 
pressure  field  and  the  associated  moments  for  liquid-filled  cylinders  embedded  in  spinning  and  nutating  j 
projectiles.  The  design  goals  are  high  efficiency,  robustness,  and  the  opportunity  of  extending  the  ! 
method  to  unsteady  problems.  The  method  uses  Chebyshev-Fourier-Chebyshev  expansions  in  the  f 
radial,  azimuthal,  and  axial  direction  and  exploits  the  symmetries  of  the  problem.  We  present  solutions  ! 
for  the  steady  motion  and  compare  with  experimental  data.  We  also  evaluate  the  performance  of  our  j 
code  in  comparison  with  other  computer  codes.  .,.  ..  j 

1.  Introduction  j 

It  is  well  known  that  spin-stabilized  shells  carrying  liquid  payloads  can  suffer  a  dynamical  instabil-  • 
ity  which  results  in  an  increased  nutation  (or  yaw)  angle  and  a  simultaneous  loss  in  spin  rate.  Labora-  j 
tory  experiments,  computational  results,  and  field  tests  indicate  that  these  phenomena  arise  from  the  j 
nutation-induced  fluid  motion  in  a  certain  range  of  small  Reynolds  numbers.  Although  in  special  cases 
this  instability  has  been  overcome  by  trial  and  error,  future  design  of  reliable  projectiles  would  profit 
from  the  opportunity  to  calculate  the  liquid  moments  and  to  account  for  these  moments  in  flight  simu¬ 
lators.  The  empirical  data  base  12  is  sparse,  however,  and  computational  methods  in  use  3'  4'  5  are  • 
rather  demanding.  ’  j 

In  previous  work,0  we  conducted  a  theoretical  analysis  which  aimed  at  the  origin  of  the  viscous  j 
despin  (negative  roll)  moment.  This  analysis  showed  that  the  deviation  from  solid  body  rotation  is  j 
governed  by  a  small  parameter  <  -----  Q  s\n0/u>  involving  the  nutation  rate  fi  ,  the  nutation  angle  0,  and  | 
the  spirt  rate  u>.  A  solution  of  the  linearized  equations  was  developed  for  a  finite-length  segment  of  an  j 
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infinitely  long  cylinder,  i.  disregarding  the  end  walls  of  the  cylinder.  Velocity  field  and  the  viscous 
components  of  the  moments  were  obtained  in  closed  form.  The  velocity  field  consists  only  of  an  axial 
component  of  order  0( c)  which  is  the  dominating  feature  of  the  fluid  motion  and  produces  a  negative 
roll  moment  of  order  O  (t  *)  owing  to  Coriolis  forces  Although  these  results  arc  in  reasonable  agree- 
iiienL  with  experimental  and  computational  data,  one  may  anticipate  modifications  of  velocity  field  and 
roll  moment  if  nonlinear  terms  are  taken  into  account.  Moreover,  yaw  and  pitch  moments  contain 
essential  contributions  of  the  pressure  ''  that  originate  from  the  turning  of  the  (low  near  the  end  walls. 
The  effect  of  nonlinearity  was  studied  1  by  using  perturbation  expansions  m  c  and  was  found  to  be 
small  The  flow  in  a  finite-length  cylinder,  however,  can  only  be  captured  by  a  computational 
approach 

The  existing  computer  codes  may  serve  for  establishing  some  basic  results  but  are  too  inefficient 
and  insufficiently  verified  for  routine  applications.  Our  analytical  work  suggests  the  use  of  a  code  that 
exploits  (i)  the  near-lmcarity  of  the  governing  equations  and  (li)  the  smoothness  of  the  solution  in  the 
relevant  range  cf  Reynolds  numbers  Therefore,  we  have  pursued  a  simple  concept  that  is  open  *.o 
numerous  refinements.  We  use  Chebyshev-Fouricr-Chcbyshcv  expansions  in  r  ,  <p ,  z  ,  respectively,  and 
convert  the  lincari/cd  equations  into  a  linear  algebraic  system  for  the  expansion  coefficients.  The  solu¬ 
tion  of  this  system  (or  any  other  solution  at  neighboring  parameters)  is  used  as  initial  approximation 
fo:  iterative  improvement  by  the  modified  Newton  method.  The  feasibility  of  this  approach  has  been 
demonstrated  8  with  a  crude  .spectral  approximation  to  the  solution.  The  version  of  the  code  reported 
here  exploits  the  diametral  symmetry  of  the  flow  and  allows  for  higher  resolution  at  modest  CPU  times. 
Tins  version  can  also  be  adapted  to  a  time- accurate  analysis  of  tbe  unsteady  problem. 


2.  Governing  Equations 

We  consider  t.iie  deviation  ,  pi  from  solid-body  rotation  m  a  nutating  coordinate  system 
7  ,  y  ,  z  where  ;  is  the  cylinder  axis  and  i  is  coplanar  with  the  two  axes  of  rotation.  (For  a  detailed 
discussion  of  the  governing  equations,  see  the  article  by  Herbert  6  1  All  quantities  are  made  nondimen- 
sional  using  a  ,  w,  and  p  for  scaling  length,  time,  and  mass,  respectively.  The  solution  depends  on  four 
noiiditmmsional  parameters:  aspect  ratio  X  =  c  ja  ,  nutation  angle  0 ,  frequency  r  —  H  /w,  and  Rey¬ 
nolds  number  /?  —  f>»ja2/ii,  where  2c  is  the  length  of  the  cylinder,  p  is  the  density,  and  p  the  viscos¬ 
ity  of  the  fluid  The  aspect  ratio  enters  the  solution  only  through  the  boundary  conditions  at  the  end 
walls  of  the  cylinder  The  motion  is  subject  to  the  no-slip  and  no-penetration  conditions  at  the  cylinder 
walls.  Therefore,  the  boundary  conditions  on  the  dc.iation  velocity  arc  homogeneous 

In  cylindrical  coordinates  r  ,q,  z,  the  equations  for  the  velocity  components  =  [vr ,  v^,u,) 
and  pressure  pj  take  the  form 
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Tr  —  -  (COS#  ,  r d  —  fsin<!»  ,  r,  —  tcqsO  ,  i  —  rsin 0  .  (2) 

The  effect  of  nutation  is  comprised  in  the  (^-periodic  force  term  ~2rrr  --  2<  r  cos^  in  the  z- 
momentum  equation  (Id).  For  t  —  0,  equations  (1)  support  the  trivial  solution  vd  =  0,  pd  =  0. 
The  system  supports  the  following  symmetries: 

vr{T  i  <*>+77 ,  -  z)  —  vr  (r  ,  <t>,  s)  .  -z)  t’c(r.6,z).  (3a,b) 

v*(r,  <f>+T7,-z)  —  -v.(r.d.s),  pd  (r  ,  6  )  pd{r.o.z).  (3c, (I) 


3.  Spectral  Approximations  for  a  Finite-Length  Cylinder 

The  results  of  the  analytical  work  suggest  that  a  good  approximation  to  the  How  in  a  Unite 
cylinder  can  be  obtained  by  solving  linearized  versions  of  equations  (l).  Linearization  can  be  performed 
in  different  ways.  The  first  is  a  linearization  in  f .  Beyond  cqs.  (3).  the  resulting  equations  support  the 
additional  symmetries 

vrf(r  ,  <t>+n  ,  2  )  =  -  vrf(r  ,  6.  z  )  ,  pd(r  ,  .  z)  -  pd(r  .  6,  :  )  .  (*l) 

These  relations  permit  useful  checks  on  the  results  of  the  spectral  code.  A  second  linear  system  can  be 
obtained  by  linearization  in  the  components  of  vd  .  This  linearization  retains  coupling  terms  such  as 
2 t i,vz  in  eq.  (lb)  which  destroy  the  symmetries  (-1).  The  second  system  can  be  considered  a  special  case 
with  vj0*  =  p^0)  —  0  of  a  third  linearization  about  some  known  solution  vj°\  This  third  pro¬ 

cedure  is  very  efficient  if  the  solution  is  sought  for  a  densely  spaced  sequence  of  parameter  combinations 
as  in  flight  simulations.  ; 

The  algebraic  form  of  the  equations  is  obtained  by  use  of  spectral  collocation.  The  velocity  com¬ 
ponents  are  expressed  in  the  form 

=  E  E  E  «Wm  Rkt(r)  F,(t)  Z,“,(r  /X)  ,  (5) 

1  =1  1  =1  m  =1 

with  similar  expressions  for  v  d>1  v, ,  and  pd  .  The  azimuthal  functions  are  F,  —  cos  i(/  -  l)<£/2]  for  odd 
Fi  —  sin  V  4> /-]  for  even  l,  where  /  =  1,  2,  •  L,  and  L  is  odd.  The  azimuthal  collocation 
points  are  equidistant,  <f>l  =  2t t{1  -  l)/L  .  The  expansion  functions  in  radial  and  axial  direction 
depend  on  the  index  /  and  may  lie  different  for  the  variables  vr  ,  v.  .  and  pd  .  They  are  combina¬ 
tions  of  even  or  odd  Chebyshev  polynomials  such  that  (i)  the  homogeneous  boundary  conditions  are 
implicitly  satisfied,  (ii)  the  symmetry  conditions  (3)  are  satisfied,  and  (iii)  the  limit  values  of  the  vari¬ 
ables  for  r  — ♦  0  (i.  e.  the  values  on  the  axis)  arc  independent  of  <$> .  The  collocation  points  are  rt  — 
cos  [(2 k  -  l)ir /4K  j,  A:  =  1,2,  •  •  •  K ,  and  zm  /X  =  cos  [(2m  -  1)77 /-LUI.  m  —  1,2,  •  •  •  M ,  Con¬ 
sequently,  0  <  rt ,  and  no  points  are  located  on  the  axis.  Also,  rk  <  1.  zm  <  X.  such  that  no  points 
are  located  on  the  surface.  This  choice  prevents  the  occurrence  of  spurious  pressure  terms  and  avoids 
the  difficulties  associated  with  corners  and  axis.  The  points  are  concentrated  near  the  boundary  such 
that  higher  resolution  in  this  region  is  obtained  without  additional  coordinate  stretching. 

The  spectral  collocation  method  converts  the  linear  system  of  partial  differential  equations  derived 
from  eqs.  (l)  into  an  algebraic  system  of  dimension  /V  —  4KLM  for  the  coefficients  vktm  ,  , 

wtim  ,  and  pklm  of  vr ,  v tj>,  v,  .  and  pd ,  respectively.  The  linear  algebraic  system  for  the  expansion 
coefficients  is  solved  by  Gauss  elimination  with  partial  pivoting.  The  subroutine  used  retains  all  data 
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required  to  solve  the  same  system  with  a  new  right-hand  side  without  repeating  the  costly  reduction  of 
the  matrix  to  upper  triangular  form.  Once  the  solution  is  obtained,  a  new  right-hand  side  is  formed 
taking  the  nonlinear  terms  into  account  and  the  system  b  iteratively  solved  until  sufficient  accuracy  is 
achieved  The  procedure  is  equivalent  to  the  modified  Newton  iteration  (without  updating  the  Jaco¬ 
bian  in  every  step)  and  converges  rapidly  since  the  nonlinear  corrections  to  the  velocity  are  small  while 
the  pressure  appears  linear  in  equations  (1). 

4.  Results  of  the  Spectral  Code 

In  the  following,  we  present  some  results  for  R  --  11.95,  0  --  -20".  r  —  0.1667.  and  X  =  4.368 
which  results  in  e  —  0.057.  The  results  are  for  /\  -  f.  -  M  5.  and  consequently  N  —  500.  Fig¬ 

ure  l  shows  the  axial  and  radial  velocities  in  the  planes  6  ~  15*  and  135*  .  Only  the  upper  half, 
z  >  0,  of  the  cylinder  is  shown;  the  lower  half  ts  governed  by  the  symmetries  (3)  The  velocity  distri¬ 
bution  at  2  --  0  agrees  well  with  the  results  of  the  perturbation  analysts  and  computations  with  the 
Sandia  code.  Near  the  end  walls,  the  solution  seems  to  be  more  realistic  and  more  accurate  than  the 
Sandia  results.  Good  resolution  of  the  velocity  gradients  and  pressure  at  the  boundary  is  important  for 
accurate  calculations  of  the  moments  The  figure  also  verifies  the  existence  of  a  predominantly  axial 
flow  over  most  of  the  cylinder  length,  except  within  a  region  of  the  order  of  the  radius  near  the  end 
wall.  Linear  and  nonlinear  velocity  distributions  are  hardly  distinguishable.  Clearly  visible  is  the  turn¬ 
ing  of  the  Row  near  the  end  wall. 

The  pressure  distribution  in  the  plane  -  45*  is  shown  in  figure  2  with  the  heavy  lines  marking 
positive  values.  Remarkable  are  the  regions  of  high  and  low  pressure  in  the  corners  near  6  45*  and 

<j>  ss  225°  ,  respectively,  which  produce  large  pressure  contributions  to  the  moments  about  x-axis  and 

y  -axis. 

I  he  dominant  components  of  velocity  and  pressure  fields  are  azimuthallv  periodic  with  period  2v 
The  harmonics  arc  small,  indicating  the  small  effect  of  nonlinearity.  The  only  important  nonlinear  term 
is  the  aperiodic  mean  flow.  The  axial  and  radial  mean  velocity  field  is  given  in  figure  3.  Thu  streaming 
term  produces  a  toroidal  mean  motion  near  the  end  wall. 


The  despin  moment  about  the  r-axis  is  largely  governed  by  the  shear  stress  at  the  side  wall 
caused  by  the  azimuthal  mean  velocity.  Thu  component  shows  little  variation  over  about  90%  of  the 
cylinder  length.  This  result  explains  the  good  agreement  of  our  earlier  analytical  results  for  the  despin 
moment  e  with  experimental  and  computational  data.  Besides  the  calculation  of  the  moments  froiti 
local  velocity  gradients  and  pressures  at  the  surface,  we  have  also  obtained  these  moments  from  volume 
integrals  involving  the  velocity  only.  This  second  method  provides  accurate  values  of  the  moments 
from  low-resolution  spectral  approximations  that  would  be  insufficient  when  using  the  surface  values. 
Reasonable  accuracy  of  the  moments  can  be  obtained  with  K  —  M  =  4.  L  —  3.  Thu  approximation 
has  been  used  to  obtain  the  data  shown  in  figure  4  together  with  results  of  the  Sandia  code  4  and  the 
analytical  results.  6  The  region  of  maximum  despin  moment  will  be  subject  to  further  study  with  higher 
resolution 


5.  Discussion 

I'he  finite-difference  code  developed  at  Sandia  Laboratories  3  **  provides  the  steady  solution  at 
11  X  24  X21  grid  points  in  r  .6.  2 -direction  by  integrating  over  typically  104  to  8  I04  time  steps,  a 
task  that  requires  6  to  48  minutes  CPU  time  on  a  Crav-lS.  This  requirement  translates  into  6  to  48 
days  on  the  MC6S010  based  Apollo  DN300  work  station  used  for  our  studies.  The  result  consists  of 
over  22,000  values  for  velocities  and  pressure.  Strik werda  Sz  Nagel  1  briefly  describe  a  code  using  finite 
differences  in  radial  and  axial  directions  and  pscmlospcctrn!  differencing  in  the  azimuthal  direction. 
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Nonuniform  grids  arc  introduced  for  increased  resolution  near  the  walls.  The  difference  equations  are 
solved  by  an  iterative  method  based  on  successive  over- relaxation.  The  computer  time  required  is  com¬ 
parable  to  that  of  the  Sandia  code  (Nusca,  BRL,  personal  communication).  The  relative  merits  of  the 
two  codes,  espec  ,  !ly  with  respect  to  the  captured  range  of  Reynolds  numbers  are  yet  concealed. 

I  lie  experience  with  the  present  version  of  the  spectral  code  shows  that  high  performance  can  be 
achieved.  The  reported  run  with  >V  —  500  requires  about  3  hours  on  an  Apollo  DN300  work  station. 
The  solution  is  obtained  in  semi-analytical  form  with  only  A'  500  numerical  coefficients.  This  low 
data  volume  is  very  attractive  for  communication  with  remote  .supercomputers.  The  code  is  well  suited 
for  vectorization,  since  practically  all  CPU  time  is  spent  on  constructing  anti  solving  an  algebraic  sys¬ 
tem.  "I'li c  code  demands  larger  memory  than  other  codes  because  til-bit  arithmetic  is  highly  recom¬ 
mended  for  spectral  methods  in  general,  and  the  algebraic  system  requires  A’ (A*  1)  wortls  of  storage. 

Nowadays,  this  memory  requirement  appears  acceptable  even  if  higher  resolution  is  desired.  Calcula¬ 
tion  of  the  moments  for  figure  4  (Af  —  1 0*2)  requires  less  than  12  minutes  per  point. 

Finally,  there  are  still  ways  to  improve  the  performance.  The  solution  process  can  be  split  into 
two  levels,  the  first  of  which  calculates  only  the  velocity  components  from  vorticity  equations  while  the 
pressure  is  obtained  a  posteriori  bv  solving  the  Poisson  equation.  For  calculating  the  moments  by 
volume  integrals,  only  the  first  step  is  required.  We  expect  that  the  next  version  of  our  rode  will  pro¬ 
vide  accurate  solutions  for  the  moments  in  less  than  1  minute  on  an  MCGS020/G8SS1  based  work  sta¬ 
tion  and  can  be  used  as  an  efficient  subroutine  in  flight  simulators. 
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ABSTRACT 

Artillery  shells  with  liquid  payloads  may  experience  a 
severe  flight  instability  owing  to  the  moments 
exerted  by  the  viscous  fluid  motion  in  the  cylindrical 
payload  container.  Incorporation  of  these  moments 
into  flight  simulators  as  a  routine  design  tool  requires 
a  highly  efficient  code  for  solving  the  Navier-Stokes 
equations.  We  describe  a  spectral  collocation  method 
which  is  based  on  Chebyshev-Fourier-Chebyshev 
expansions  in  the  radial,  azimuthal,  and  axial  direc¬ 
tion.  The  method  exploits  the  symmetries  of  the 
problem.  Using  a  volume  approach  and  an  analytical 
result  by  Rosenblat,  accurate  moments  are  obtained 
in  small  fractions  of  the  time  required  by  other 
codes.  Solutions  for  the  steady  motion  are  presented 
and  compared  with  numerical  and  experimental  data. 
The  performance  of  our  code  is  evaluated  in  com¬ 
parison  with  other  computer  codes. 

Introduction 

Gyros  and  rotating  fluids  often  exhibit  unex¬ 
pected  behavior.  In  the  past,  it  has  been  recognized 
that  spin-stabilized  shells  with  liquid  payloads  can 
suffer  a  dynamical  instability  originating  from  reso¬ 
nance  with  inertial  waves.1  Since  this  phenomenon  is 
basically  inviscid  and  is  routinely  avoided  by  proper 
design,  it  was  surprising  to  observe  in  some  cases 
another  type  of  instability  which  is  characterized  by 
an  increase  in  nutation  (or  yaw)  angle  and  a  simul¬ 
taneous  loss  in  spin  rate.  The  rapid  drop  in  spin  rate 
is  clearly  a  viscous  phenomenon,  and  laboratory 
experiments,  computational  results,  and  field  tests 
have  meanwhile  shown  that  this  instability  is  caused 
by  the  nutation-induced  fluid  motion  in  a  certain 
range  of  relatively  small  Reynolds  numbers. 
Although  in  special  cases  this  instability  has  been 
overcome  by  trial  and  error,  future  design  of  reliable 
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projectiles  would  take  profit  from  the  opportunity  to 
calculate  the  liquid  moments  and  to  account  for 
these  moments  in  flight  simulators.  The  empirical 
data  base  ?  is  sparse,  however,  and  computational 
methods  in  use  3- 4’ 5' 6  are  rather  demanding.  An 
evaluation  and  verification  of  the  codes  by  Vaughn  et 
at.  4  and  Strikwerda  &  Nagel  5  is  currently  conducted 
at  BRL.7  Typical  computer  times  for  a  single  case 
are  in  the  range  of  6-12  hours  on  VAX-class 
machines.  Six-degree-of-freedom  flight  simulators  8 
typically  use  2  106  time  steps  over  the  flight  time  of 
the  order  of  30  seconds.  Study  of  the  interaction  of 
the  interior  fluid  motion  with  the  exterior  aeroballis- 
tics  consequently  requires  either  a  very  fast  subrou¬ 
tine  for  calculating  the  liquid  moments  or  interpola¬ 
tion  in  a  multi-dimensional  table  of  500-1000  8  pre¬ 
calculated  values.  Hence,  Sight  simulations  for 
liquid-filled  shells  are  currently  a  very  expensive  tool 
and  are  not  ready  for  routine  applications. 

In  previous  work  9  we  conducted  a  theoretical 
analysis  which  aimed  at  the  origin  of  the  viscous  des- 
pin  (negative  roll)  moment  in  cylinders  of  large 
aspect  ratio.  This  analysis  showed  that  the  deviation 
from  solid  body  rotation  is  governed  by  a  small 
parameter,  <  =  (0  /w)sin0,  involving  the  nutation 
rate  D  ,  the  nutation  angle  9,  and  the  spin  rate  w.  A 
solution  of  the  linearized  equations  was  developed  for 
a  finite-length  segment  of  an  infinitely  long  cylinder, 
i.  e.  disregarding  the  end  walls  of  the  cylinder.  Velo¬ 
city  field  and  the  viscous  components  of  the  moments 
were  obtained  in  closed  form.  The  velocity  field  con¬ 
sists  only  of  an  axial  component  of  order  O  (« )  which 
is  the  prominent  feature  of  the  fluid  motion  and  pro¬ 
duces  a  negative  roll  moment  of  order  0(<2)  owing 
to  Coriolis  forces.  Although  this  roll  moment  is  in 
reasonable  agreement  with  experimental  and  compu¬ 
tational  data,  the  original  analysis  accounted  only  for 
the  viscous  part  of  the  yaw  and  pitch  moment. 
These  latter  moments  contain  essential  contributions 
of  the  pressure  4  that  originate  from  the  turning  of 
the  flow  near  the  end  walls  and  were  not  captured  by 
the  linear  analysis  The  effect  of  nonlinearity  was 
studied  10  by  using  perturbation  expansions  in  <  and 
w  as  found  to  be  small  except  for  an  aperiodic  stream- 
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iug  term  in  the  azimuthal  direction. 

Athough  the  perturbation  approach  cannot  cap¬ 
ture  the  pressure  field,  it  provides  valuable  insight 
into  the  structure  of  both  equations  and  solution. 
The  analytical  work  suggests  the  use  of  a  numerical 
method  that  exploits  (i)  the  near-linearity  of  the 
governing  equations  and  (ii)  the  smoothness  of  the 
solution  in  the  relevant  range  of  Reynolds  numbers. 
We  have  therefore  pursued  a  simple  concept  that  is 
open  to  further  refinements.  We  use  Chebyshev- 
Fourier-Chebyshev  expansions  in  r  ,  <p ,  z  ,  respec¬ 
tively,  and  convert  the  linearized  equations  into  an 
algebraic  system  for  the  expansion  coefficients. 
Linearization  can  be  performed  about  the  trivial 
solution  or  any  other  known  solution,  e.  g.  at  neigh¬ 
boring  parameters.  The  solution  of  the  linear  alge¬ 
braic  system  is  u»ed  as  initial  approximation  for 
iterative  improvement  by  the  modified  Newton 
method.  The  feasibility  of  this  approach  has  been 
demonstrated  11  with  a  crude  spectral  approximation 
to  the  solution.  Problems  in  calculating  the  pressure 
that  arise  from  the  invalidity  of  the  basic  equations 
in  the  corners  joining  the  flat  end  walls  to  the 
cylindrical  side  wall  have  meanwhile  been  over¬ 
come.12  The  present  version  of  the  code  exploits  the 
diametral  symmetry  of  the  flow  about  the  center  of 
the  cylinder  and  allows  for  higher  resolution  at  mod¬ 
est  CPU  times.  This  version  can  also  be  adapted  for 
the  analysis  of  unsteady  problems.  Dramatic 
increase  in  efficiency  has  recently  been  achieved  13  by 
combining  an  analytical  result  of  Rosenblat  ft  al.  6 
with  a  volume  formulation  for  calculating  the  liquid 
moments.  The  moments  can  be  obtained  from  only 
the  simply  periodic  components  of  the  axial  velocity 
and  the  azimuthal  streaming  term.  A  fast  subroutine 
for  flight  simulations  exploits  the  analytical  results. 
For  more  accurate  studies,  complete  tables  of 
moments  can  be  calculated  in  a  few  hours  on  a 
VAX-type  computer. 

Governing  Equations 

We  consider  the  steady  motion  of  a  fluid  of 
density  p  and  viscosity  /z  in  a  cylinder  of  radius  a 
and  length  2c  in  an  aeroballistic  coordinate  system 
x  ,  y  ,  x  ,  where  x  is  the  axis  of  the  cylinder,  as  shown 
in  Figure  1.  The  inertial  axis  Z  in  (light  direction 
and  the  z-axis  enclose  the  nutation  angle  8.  The 
cylinder  rotates  with  the  spin  rate  w  about  z  while 
the  z,z -plane  rotates  with  the  nutation  rate  Q 
about  the  Z -axis.  Spin  rate  w  and  nutation  late  P. 
are  constant.  All  quantities  are  made  nondimen- 
sion ai  us.ng  a  ,  w,  tend  p  for  scaling  length,  tune,  and 
inass,  respectively  The  solution  depends  on  four 
nondimensional  parameters:  aspect  ratio  t)  —  c  ,  <*  , 


nutation  angle  9,  frequency  r  =»  0  /w.  and  Reynolds 
number  Re  =  pwa2//i.  The  aspect  ratio  enters  the 
solution  only  through  the  boundary  conditions  at  the 

mill  walls  of  the  cylinder.  The  intlion  is  subject  to 

the  no-slip  and  no-pcnctration  conditions  at  the 
cylinder  walls.  Since  the  velocity  field  degenerates  for 
cither  w  =  0,  Cl  —  0,  6  =  0,  or  p— »oo  to  rigid- 
body  rotation  of  the  fluid,  it  is  appropriate  to  con¬ 
centrate  on  the  deviation  v*  of  the  velocity  from 
rigid-body  rotation  vf , 

v  =  vr  +  v*  ,  v'  =  re^  ,  (1) 

where  is  the  azimuthal  unit  vector.  The  boundary 
conditions  on  v*  are  homogeneous.  The  pressure 
field  is  split  according  to 

P  =  P'  +  P *  ,  (2a) 

p'  =  ^(l+r.f  +  r^+xV-  2rz<rf  ,  (2b) 


where  r,  —  -t  cosip ,  —  <sin<£,  r,  =  rcosfl, 

(  =  rsintf.  The  pressure  pr  differs  from  the  pres¬ 
sure  in  rigid-body  rotation.  The  form  of  pr  is 
chosen  such  that  the  reduced  pressure  p *  appears 
only  in  the  z-momentum  equation. 

In  cylindrical  coordinates  r  ,  <p,  z ,  the  equations 
for  the  velocity  components  vd  —  (vr  ,  v^,  v, )  and 
pressure  p  4  take  the  form 


Id,  . 

TS1"'1 


D'v,  - 


v*2 


1  dvt  dv, 

+  7  ~dT+  TT  “ 0  ’ 

-  2(1  +  t,)v4  +  2 t4v,  = 


(3a) 


dp4 

dr 


+ 


2  dv*, 

r2  d*  1  ’ 


(3b) 


D'v.  + 


Vr  V 


r  v* 


+  2(1  4-  r ,  )vr  -  2rrt>,  = 


+  TT  +  TTir1 '  (3c) 

D’v,  +  2rrvp  -  2t4v,  = 

-^T-~  2rrr  +  -±-D"v,  ,  (3d) 

dz  Re 


where 


D ' 


D 


d  d  d_  v*  d 

dt  +  dip  "*  Vf  dr  r  dip 

•  =JL  +  1±  +  ±JL 

dr 2  r  dr  r2  dtp2 


+  v. 


d_ 
dz  ’ 


d 2 

dz2  ‘ 


The  primary  effect  of  nutation  is  the  ^-periodic  force 
term  -  2r  r,  =  2trcos<p  in  the  z-momentum  equa¬ 
tion  (3d).  For  <  =  0,  equations  (3)  have  the  trivia! 
solution  v*  «=  0,  p4  ss  0,  The  system  supports 

the  following  symmetries: 
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e,  (r  ,  <p+»r, 

-  :)  =  vr{r  ,  4>,  z)  , 

(4a) 

,  9-x. 

-  z  )  =  u*(r  ,  0,  z)  , 

(4b) 

v,  (  r  .  </>  .  ”  , 

J  )  —  ...  (  r  ,  )  , 

(■•) 

f>‘(r  , 

“  *  )  —  P*(r.  4>,  O 

(4d) 

Some  Analytical  Results 

The  steady  flow  in  a  relatively  long  cylinder 
(aspect  ratio  r/  >  4)  at  low  Reynolds  number  is 
expected  to  exhibit  little  axial  variation  over  most  of 
the  cylinder  length.  Previous  work  9  has  therefore 
relaxed  the  boundary  conditions  at  the  end  walls  and 
studied  the  steady  flow  in  a  finite  segment  of  an 
infinitely  long  cylinder. 

In  the  physical  situations  of  interest, 
c  —  (0  /w)sin0  is  a  small  parameter,  and  conse¬ 
quently,  it  is  reasonable  to  pursue  a  perturbation 
expansion  in  e  .  This  provides  in  the  form 

v'  -  (5) 


and  similar  expressions  for  p  d  .  The  development  of 
expressions  for  the  expansion  coefficients  v*n  *  from 
equations  (3)  leads  to  an  alternating  pattern 

I  (0.  0,  ty"!)  ,  n  odd  , 

V*  ^  I  (ur[,>,  v4  (">,0)  ,  n  even  . 
and  the  components  of  v**)  take  the  form 

+  ume-iim4)  ,  (Ta) 


vi"1  =  v„  +  £  .  (7b) 

m  —  1 

Vj(n)  +  t2>lm  e -.(2«- II#)  .(7c) 


where  the  tilde  denotes  the  complex  conjugate  The 
aperiodic  term  in  v f,"i  is  suppressed  by  the  con¬ 
tinuity  equation  The  r-dependent  coefficient  func¬ 
tions  in  Eqs.  (7)  are  required  to  satisfy  homogeneous 
boundary  conditions  at  r  =  1  and  to  be  finite  at  the 
axis  r  —  0. 

The  axial  velocity  at  order  0(<)  can  be  found  in 
analytical  form, 

/,/ar  ) 

«":i(r)  —  '  !~ ,  -  r  i  •  W 

i  i(f-») 


where  /)  denotes  the  modified  Bessel  function,  and 
a  =  (1  -»•  i)(fie  2)1, 7  This  solution  b  valid  for 
arbitrary  Reynolds  number  but  may  be  unstable  as 
H(  exceeds  some  critical  value.  This  component  is 
the  dominating  feature  of  the  flow  in  a  long  cylinder 
The  interesting  properties  of  the  associated  flow  field 


are  discussed  by  Herbert.9 

At  order  0(e2),  comparison  of  the  equation  for 
v;n  with  the  imaginary  part  of  the  equation  for  ui  n 
immediately  shows  that  the  aperiodic  component  oT 

the  azimuthal  velocity  is 

t'ao(f)  =  -2  Iin(to  i  ,(r  )|  .  (fi) 

'flic  0-pcriodic  components  are  governed  by  a  cou¬ 
pled  set  of  inhomogeneous  differential  equations  with 
variable  coefficients  Essential  simplification  at  the 
expense  of  increasing  the  order  of  differentiation 
results  from  eliminating  by  use  of  the  continuity 
equation.  With  some  effort,  the  radial  velocity  com¬ 
ponent  of  0  (<2)  can  be  found  in  closed  f  'r:n. 

“2i(*  )  —  )| 

««  2i  v2  J 'A*  />/?) 

3  a3  »  TJiTfJZ) 

where  a  =  0r  ,  /?  =  (i  -  1  )Rtx^,  and  J J  2,  and 
Y-2  arc  Bessel  functions.  The  coefficients  cj,  c  2.  tj, 
and  c,  can  be  determined  numerically.10  In  view  of 
the  effort  involved  in  deriving  the  closed  form  solu¬ 
tion  for  u21  and  the  ultimate  need  to  determine  the 
coefficients  numerically,  the  differential  equations  for 
the  third-order  components  were  solved  by  means  of 
a  spectral  collocation  method. 

The  motion  is  governed  by  the  axial  component 
w  1  j  at  order  0(c).  Of  the  higher  order  terms,  only 
the  aperiodic  term  t>20  13  substantial.  In  the 

cylinder’s  center  section,  these  terms  are  in  good 
agreement  with  computational  results.  Ail  the  other 
terms  are  not  only  of  order  0(1)  but  in  fact  less  than 
unity,  assuring  rapid  convergence  of  the  perturbation 
series  The  contribution  of  uijj  to  the  despin 
moment  is  negligible  The  d>-periodic  terms  oscillate 
about  zero  as  r  varies  between  0  <  r  <1.  Accu¬ 
rate  representation  of  single  high-order  terms  by 
radial  Chebyshev  series  may  require  numerous  expan¬ 
sion  functions  For  the  total  velocity  field,  however, 
the  error  in  representing  these  terms  is  of  little 
importance.  At  Reynolds  numbers  in  the  range  of 
maximum  despin  moment,  reasonably  accurate 
approximations  can  be  obtained  with  only  a  few 
polynomials  in  radial  direction.  In  the  azimuthal 
direction,  the  solution  is  governed  by  terms  periodic 
m  o,  and  by  the  aperiodic  term  Fourier  series 

with  three  or  five  modes,  therefore,  provide  approxi¬ 
mations  of  sufficient  accuracy  for  practical  purpose. 
The  perturbation  analysis  clearly  shows  that  the 
main  features  of  the  flow  are  governed  by  the  linear 
O  (<  I  part  of  equations  (3)  w  ith  small  corrections  for 
nonlinearity  This  property  will  not  change  for  a 
finite-length  cylinder 
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Spectral  Approximations 

The  results  of  the  analytical  work  suggest  that 
a  good  approximation  to  the  flow  in  a  finite  cylinder 
can  be  obtained  by  solving  linearized  versions  of 
equations  (3).  l,ip«anz*tinn  ran  he  performed  in 
different  ways.  The  first  is  a  linearization  in  < 
Besides  Eqs.  (4),  the  resulting  equations  support 
additional  symmetries: 

v'(r,6+!r,;)-"-v'(r1«.:)1  (Mi) 

p<(r,<i4-!r.:)  =  -p<(r1ip,:).  (Mb) 

These  relations  permit  useful  checks  on  the  results  of 
the  spectral  code.  A  second  linear  system  can  be 
obtained  by  linearization  in  the  components  of  v* 
This  linearization  retains  coupling  terms  such  as 
2 r4v,  in  Eq  (3b)  which  destroy  the  symmetries  (11). 
The  second  system  can  be  considered  a  special  case 
with  =  0  of  a  linearization  about  aome  known 
solution  v'  The  latter  procedure  is  very  efficient  if 
tb“  solution  is  sought  for  a  densely  spaced  sequence 
of  parameter  combinations  r~  in  flight  simulations 

The  algebraic  form  of  the  equations  is  obtained 
by  use  of  a  spectral  collocation  method.  The  velocity 
components  are  expressed  in  the  form 

i',  =  E  £  £  nkltnRi(r)Fl(o)Zm(^)  ,  (12) 

*  -  1  I  =  1  m  -  |  'I 

with  similar  expressions  for  r#,  t2 ,  and  p i  The 
azimuthal  functions  are 


F, 


cos  - <?  ,  /  odd  , 

2 

/  j.  , 

sin  ~i>  ,  l  even  . 


(13) 


The  azimuthal  collocation  points  are  equidistant. 


=  2 -(/  -  1  ),L  ,  /  =  1,2.  •  L  .  (14) 


and  L  is  odd. 

In  a  first  version  of  the  code,  radial  and  axial 
collocation  points  are  located  at  the  maxima  of  the 
highest  Chebyshev  polynomials  The  boundary  con¬ 
ditions  are  implemented  by  replacing  three  of  the 
four  differential  equations  in  the  boundary  points. 
The  question  then  is  which  equation  should  be 
retained  and  where  the  condition  on  the  pressure, 
e  g  p  d  =  0.  should  be  applied.  Trial-and-error 
leads  to  numerous  cases  with  ill-determir  :d  matrices 
or  zero  determinant.  In  other  cases,  a  correct  solu¬ 
tion  for  the  velocity  field  is  obtained,  but  the  pres¬ 
sure  contains  a  non-physical  spurious  term  Prob¬ 
lems  with  specli.il  calculations  of  the  pressure  in 
closed  domains  with  corners  are  well-known  hut  the 
reports  on  their  origin  and  methods  for  solution  are 


rather  unspecific.  We  have  therefore  performed  a 
detailed  analysis  of  the  flow  in  a  square  driven  by  an 
internal  force  field.  This  simpler  two-dimensional 

1 1 ruble i ii  exhibits  all  cliuiac  tcri«tic«  -  including  lltr 

spurious  pressure  term  -  of  the  original  problem. 
Detailed  results  of  this  study  will  be  reported  else¬ 
where.  **  The  study  reveals  that  the  spurious  term 
vanishes  in  all  collocation  points  except  the  corners, 
where  it  may  assume  arbitrary  values.  The  term  can 
be  suppressed  by  retaining  in  the  corners  one  of  the 
momentum  equations  that  contain  the  derivative  of 
the  pressure  in  the  direction  of  the  boundary. 

In  a  second  version  of  the  spectral  code,  the 
problems  of  the  pressure  calculation  have  been 
avoided  by  using  a  different  set  of  collocation  points. 
The  expansion  functions  in  radial  and  axial  direction 
depend  on  the  index  /  and  may  be  different  for  the 
variables  vr ,  v  4,  v, ,  and  p  1 .  They  are  combinations 
of  even  or  odd  Chebyshev  polynomials  such  that 

(i)  the  homogeneous  boundary  conditions  arc 
implicitly  satisfied, 

(ii)  the  symmetry  conditions  (4)  are  satisfied,  and 

(in)  the  limit  value  of  the  variables  for  r  — *  0  (i.  e. 
the  value  on  the  axis)  is  indepenH»nt  of 

The  collocation  points  are 

f*  =  S'n  k 2K  n  ’  k  ~  l'2,  K  ,  (15a) 

— —  =-  sin  —  -  - -  -  n  ,  m  =1,2,  ■  -  ■  M  .  (15b) 
r;  2  Af 

Consequently  0  <  r4  ,  and  no  points  are  located  on 
the  axis.  Also,  r*  <  1,  zm  <  r?  such  that  no  points 
are  located  on  the  surface.  The  points  in  radial  and 
axial  direction  are  concentrated  near  the  boundary 
such  that  high  resolution  in  this  region  is  obtained 
without  additional  coordinate  stretching.  Thus  the 
boundary  layers  forming  at  higher  Reynolds  number 
can  be  resoived  by  slightly  increasing  K  and  M . 

The  spectral  collocation  method  converts  the 
linear  system  of  partial  differential  equations  derived 
from  Eqs.  (3)  into  an  algebraic  system  of  dimension 
.V  =  I  K  I  M  for  the  coefficients  uWm  ,  vilm  ,  wilm  , 
and  pktm  of  vr  ,  v4,  v,  ,  and  pi,  respectively.  The 
linear  system  for  the  expansion  coefficients  is  solved 
by  Gauss  elimination  with  partial  pivoting,  lhe  sub¬ 
routine  used  retains  all  data  required  to  solve  the 
same  system  with  a  new  right-hand  side  without 
repeating  the  costly  reduction  of  the  matrix  to  upper 
triangular  form.  Once  the  solution  is  obtained,  a 
new-  right-hand  side  is  formed  by  taking  the  non¬ 
linear  terms  into  account  and  the  system  is  itera¬ 
tively  solved  until  sufficient  accuracy  is  achieved 
The  procedure  is  equivalent  to  the  modified  Newton 
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iteration  (without  updating  the  Jacobian  in  every 
step)  and  converges  rapidly  since  the  nonlinear 
corrections  to  the  velocity  are  small  while  the  pres¬ 
sure  appears  linear  in  equation*  (11) 

Resulto  for  Velocity  and  Pressure 

In  the  following,  wo  present  some  rrsulls  for 
the  velocity  and  pressure  fields  at  0  =  ‘JO1’  , 
t  ~  0.16667,  and  r)  —  4.368  which  results  in 
t  —  0.057.  The  results  are  for  K  =  6,  L  —  5,  and 
M  =  8,  and  consequently  ,Y  —  960.  Calculation  of 
a  single  solution  with  this  high  resolution  requires 
about  2  minutes  on  a  Cray-lS.  Figure  2  shows  the 
axial  and  radial  velocities  in  the  planes  o  =  45°  and 
<p  =  135°  at  Re  =  20  Only  the  upper  half. 
z  >0,  of  the  cylinder  is  shown;  the  lower  half  is 
governed  by  the  symmetries  (4).  The  scale  values 
give  the  velocity  per  unit  length  where  the  diameter 
is  six  units.  The  velocity  distribution  at  z  —  0 
agrees  well  with  the  results  of  the  perturbation 
analysis  and  computations  with  the  Sandia  code.5 
Near  the  end  walls,  the  solution  is  more  realistic  and 
more  accurate  than  the  Sandia  results  The  figure 
also  verifies  the  existence  of  a  predominantly  axial 
flow  over  most  of  the  cylinder  length,  except  within  a 
region  of  die  order  of  the  rauius  r.ea*  the  end  wall 
Linear  and  nonlinear  velocity  distributions  are  hardly 
distinguishable.  Clearly  visible  is  the  turning  of  the 
(low  near  the  end  wall.  While  the  flow  appears 
steady  in  the  coordinate  system  chosen,  the  velocity- 
field  describes  in  fact  an  oscillatory  motion  of  fluid 
elements  about  their  near-circular  orbit. 

The  pressure  distribution:  for  the  same  case  are 
shown  in  Figure  3  with  the  heavy  lines  indicating 
positive  values.  Remarkable  is  the  formation  of 
regions  of  high  and  low  pressure  in  the  corner  near 
Q  ~  45°  and  ■?  **  135°  respectively,  which  produce 
large  contributions  to  the  moments  about  z  axis  and 
y  -axis.  Except  in  this  region  near  the  end  walls,  the 
variation  of  the  pressure  is  relatively  weak  The 
azimuthal  position  of  extremum  pressure  changes 
from  o  —  0  for  small  values  of  He  to  q  ~  90°  at 
He  =  HiOO 

The  dominant  components  uf  velocity  ami  pres¬ 
sure  fields  are  azimuthally  periodic  with  period  2rr 
The  liarmoiii'-s  are  small,  indicating  the  smail  effect 
of  nonlinearity  in  the  range  of  low  Reynolds 
numbers.  The  or!,,  important  nonlinear  term  is  the 
aperiodic  mean  flow  This  is  clearly  shown  by  Figure 
4  width  gives  the  azimuthal  velocity  in  the  renter 
plane  z  0.  The  aperiodic  component  is  opposite 
to  t!ii-  rigid-body  rotation  and  exerts  a  negative  roil 
moment  li.iuugh  the  wail  shear  stress  rri  The  axial 


and  radial  mean  velocity  field  is  given  in  Figure  5. 
This  streaming  term  exhibits  a  toroidal  motion  near 
the  end  in  each  half  of  the  cylinder  and  causes  a  slow 
drift  uf  lluiil  rlriiu’iitn  with  rnaport  to  circular  orbit* 
This  mean  velocity  produces  the  symmetric  pattern 
in  How  visualizations  18  at  low  Reynolds  numbers. 

At  tin-  higher  Reynolds  number  Re  —  300,  the 
maximum  axial  velocity  appears  at  d>  ~  90°  .  /Vs 
shown  in  Figure  6,  the  flow  in  the  plane  4>  =  90° 
breaks  up  into  two  swirls,  one  in  each  half  of  the 
cylinder,  with  little  flow  across  the  plane  z  =  0. 
Three  weak  swirls  develop  in  the  plane  4>  —  0  such 
that  the  velocity  field  is  reminiscent  of  a  chain  with 
five  links.  Notably,  the  break-up  into  cells  is  res¬ 
tricted  to  an  inner  region  of  the  cylinder.  The 
motion  in  the  pronounced  boundary  layer  visible  in 
the  plane  d>  =  0  does  not  follow  the  cellular  struc¬ 
ture  and  may  have  a  direction  opposite  to  the  coie 
flow.  The  pressure  variation  is  characteristically 
different  from  that  at  low  Reynolds  number.  Figure 
7  shows  the  strong  variation  and  the  formation  of  an 
almost  symmetric  pattern  along  the  cylinder  in  the 
plane  o  =  0,  while  the  variation  at  <p  =  90°  is 
rather  weak.  This  pressure  field  explains  the  void 
observations  of  Miller  18  which  show  a  wavy  distor¬ 
tion  of  the  void  in  the  plane  <t>  —  0  at  high  Reynolds 
numbers  The  steep  and  opposite  pressure  gradients 
across  the  cylinder  axis  near  z  /rj  =  0  25  and 
z  r 7  =  0  75  displace  the  void  near  these  positions  in 
opposite  directions  along  the  diameter  at  <t>  ~  -  15° 

Calculation  of  the  Liquid  Moments 

Conservation  of  angular  momentum  for  the 
steady  flow-  in  a  control  volume  V  with  surface  5 
rotating  with  constant  rate  fi  about  a  fixed  axis 
requires 

M  ~  /  (r  X  F)d5  =  /  r  X  (20  xv)pdV 
s  v 

+  /  r  X  ;n  X  (fl  X  r)! pdV 
v 

~  f  (r  X  v)p(vd  S)  (16) 

s 

where  the  velocity  v  is  measured  relative  to  the  aero- 
b.illislR  frame  On  the  left-hand  side,  M  is  the  resul¬ 
tant  torque  on  the  control  volume,  r  is  the  position 
vector,  and  F  the  stress  acting  on  the  cylinder.  The 
presence  and  meaning  of  certain  terms  depend  on  the 
choice  of  the  control  volume.  The  surface  integral  on 
tfie  right-hand  side  of  Eq.  (16)  vanishes  if  the  surface 
of  the  control  volume  is  closed. 

Foi  ease  of  practical  application,  we  express  the 
moment  M  =  (Af,  .  A ft.\f.)  in  terms  of  cartesian 
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•  ■■  m-s  wincii  piovide  yaw,  pitch,  and  roll 

moment,  respectively.  Analogue  to  Eq  (1)  we 
decompose  the  moments  into 

M  =  M'  +  M*  ,  (t?) 

where  M'  corresponds  to  the  pure  rigid-body  motion 
while  M*1  originates  from  the  deviation  velocity  .-uni 
pressure  For  the  cylindrical  control  volume,  tin- 
rigid-body  rotation  causes  only  a  pitch  component 

,v;-s»«,n  +  ^r(i-  |V)|.  (.8) 

while  ,\/,r  -=  M,r  =  0.  Note  that  M  is  dimensionless, 
the  reference  moment  is  pasw2. 

The  evaluation  of  the  components  of  M*  bears 
some  ambiguity  that  can  be  exploited  for  advantages. 
Previous  computational  work  ^  6  employed  a  con¬ 

trol  volume  consisting  of  an  “empty”  closed  cylinder 
with  only  the  pressure  and  stresses  acting  on  me 
inside  surface  In  this  “surface  approach,”  the  right 
hand  side  of  Eq  (16)  vanishes  and  the  moments  are 
obtained  from  the  stresses  F  at  the  inside  wall  of  the 
control  volume.  Here,  we  choose  a  “volume 
approach”  that  bears  great  advantages  especially  for 
computational  work. 

We  consider  a  control  volume  consisting  of  a 
solid  cylindrical  surface  completely  enclosing  the 
liquid.  The  moment  calculation  for  this  “full"  con¬ 
trol  volume  rests  on  the  relation 

M'  -  J  r  X  (20  Xv*)pdV  (10) 

v 

Using  analytical  relations  derived  by  Rosenblat  et 
al.,6  the  components  of  M*  can  be  shown  to  take  the 
form 


St*  =■■  -  /,  ccs0 

(20a) 

Sty  —  /2  sin#  -  / 3  cost) 

(20b) 

St*  /«  sintf 

(20c) 

j  r(urcos<?  -  v  #sin<?  )rdri  <t>  dz  , 

V 

(21a) 

=  —■  f  v  p  r  ~drd  <?  dz  , 

-  V 

(21b) 

/ 3  =  -  j  !■_.  sin (ir7drdodz  , 

V 

(21c) 

/,  =---/,. 

(2  Id) 

Finally,  we  obtain  the  moments  in  the  form 

r  if  1 

St*  =  -'—n  f  J  f  v,r2cos4,Jrd<?dz  ,  (22a) 

-  .j  0  0 

*3  1 

A/„’  -  (  f  j  j  v  sr-  Jr d  v  dz 

"  i  0 


r)  1  2'  1 

- - — —  J  f  f  v,  r2sm<p  drd  4>  di  ,  (22b) 

lanS  -  ft  o  ‘6 

K(f  —  Af/lan 0  .  (22r) 

The  volume  integral  approach  leads  to  handy  expres¬ 
sions  which  involve  only  tlie  radial  and  azimuthal 
velocity  components.  Integration  over  4>  reduces  the 
requirements  in  fact  tr  the  knowledge  of  the 
aperiodic  component  of  u  *  and  the  simply  periodic 
components  of  tq  .  Therefore,  the  volume  approach 
can  also  be  applied  to  the  analytical  results  given 
above  and  provides  yaw  and  pitch  moments  without 
explicit  knowledge  of  the  pressure. 

Results  for  the  Liquid  Moments 

While  velocity  and  pressure  fields  are  primarily 
of  basic  fluid  mechanical  interest,  the  practical  need 
for  the  moments  dictates  the  measure  for  efficiency 
of  the  code.  The  moments  derived  from  the  volume 
approach  and  the  surface  approach  applied  to  the 
same  spectral  solutions  is  shown  in  Tables  1  and  2, 
respectively.  The  Reynolds  number  Re  =  20  is  in 
the  range  of  maximum  despin  moment  St,  . 


Table  1.  Volume  Approac 

h 

II 

c-  2c 

4.368  t 
L  St 

=  .1667 
St, 

e  =  20 

...  M, 

Rc  =  20 
St, 

3 

3 

3 

0.08305 

0.07475 

0.03023 

4 

3 

4 

0.08260 

007334 

0  03006 

5 

3 

5 

0.08300 

0.07332 

0.03021 

5 

3 

6 

0.08317 

0.07353 

0.03027 

6 

3 

5 

0.08300 

0.07332 

0.03021 

6 

3 

6 

0  08317 

0.07353 

0.03027 

4 

5 

4 

0.08280 

0  07353 

0.03014 

5 

5 

5 

0.08322 

0  07355 

0  03029 

6 

5 

6 

0.08340 

0.07374 

0  03035 

6 

5 

8 

0.08335 

0.07385 

003034 

It  is  obvious  that  the  volume  approach  provides 
results  of  superior  quality  and  more  rapid  conver¬ 
gence  The  required  (absolute)  accuracy  of  10' 3  for 
engineering  applications  can  be  achieved  with  the 
low  truncation  K  —  4,  L  —  3,  St  =  4.  This  accu¬ 
racy  has  to  be  seen  in  the  light  of  considerable  uncer¬ 
tainty  in  the  moments  governing  the  exterior  aero¬ 
dynamics  of  the  projectile.  We  note,  however,  that 
an  increase  in  the  aspect  ratio  may  require  additional 
expansion  functions  in  axial  direction  while  increas¬ 
ing  Reynolds  number  require:-  higher  resolution  in 
both  radial  and  axial  direction  f  igure  9  gives  the 
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Table  2,  Surfac 

e  Approach 

n 

=  4.368 

.1667 

0  —  20  Ro  -  20 

K 

L 

M 

M, 

My 

M. 

3 

3 

3 

0.07394 

0.09396 

0.03308 

4 

3 

4 

0.07247 

0.08133 

0.0201)2 

5 

3 

5 

0.07904 

0.07291 

0.03024 

5 

3 

6 

0.08178 

0.07039 

0.03028 

6 

3 

5 

0.07864 

0.07354 

0.03023 

6 

3 

6 

0.08137 

0.07115 

0.03027 

4 

5 

4 

0.07289 

0.08354 

0.02999 

5 

5 

5 

0.07894 

0.077CC 

0.03032 

6 

5 

6 

0.08152 

0.07491 

0.03036 

6 

5 

8 

0.08289 

0.07415 

0.03034 

comparison  of  the  roll  moments  for  a  wide  range  of 
Reynolds  numbers  with  the  experimental  results  of 
Miller  2  and  with  computational  results.4'0  The  devi¬ 
ation  of  the  results  of  the  Sandia  code4  is  due  to 
using  inappropriate  formulas  for  the  moments  in  the 
nutating  coordinate  system.7  The  agreement  with  the 
other  computational  data  is  good.  Test  runs  with 
high  resolution  up  to  Re  —  300  indicate  that  the 
small  errors  are  due  to  insufficient  resolution  of  the 
finite-element  code  0  in  combination  with  the  surface 
approach  for  the  moments.  The  experiments  were 
made  in  a  range  of  spin  rates  u  between  2000  and 
4000  rpm.  While  w  =  3000  rpm  has  been  used  in 
figure  9,  assumption  of  a  lower  value  would  improve 
the  comparison  with  respect  to  the  maximum  values. 
Figure  10  shows  a  similar  comparison  for  the  yaw 
and  pitch  moments.  The  results  of  the  Sandia  code 
arc  suppressed  since  they  suffer  from  a  dimensional 
inconsistency.7  While  the  agreement  for  the  yaw 
moment  at  high  Reynolds  numbers  is  surprisingly 
good,  the  deviation  in  the  pitch  moment  is  likely  to 
originate  from  insufficient  resolution  of  the  steep 
pressure  gradients.  This  source  of  discretization 
errors  has  been  eliminated  in  the  results  of  the  spec¬ 
tral  code. 

The  effect  of  nonlinearity  in  c  on  the  yaw  and 
pitch  moments  is  shown  in  Figure  11  for  different 
nutation  angles  0  at  fixed  r  =  1  and  different  values 
of  r  at  0  —  20°  .  The  figure  shows  the  ratio  of  the 
nonlinear  to  the  linear  solution.  In  absence  of  non¬ 
linearity  all  points  should  be  located  on  the  horizon¬ 
tal  lines.  Yaw  and  roll  moments  are  largely  indepen¬ 
dent  of  c.  However,  the  nonlinear  effect  on  the  pitch 
moment  is  only  negligible  at  small  nutation  angles  0 
as  they  may  occur  for  soft  launch  conditions  and 


stable  projectile  flight. 

Figure  12  shows  the  dependence  of  the  yaw  and 
pitch  moments  per  unit  length  (the  roll  moment  is 
proportional  to  A-/,)  on  the  aspect  ratio  of  Lite 
cylinder  and  compares  witli  results  of  the  code  writ¬ 
ten  hy  Strikwenla  &  Nagel.0'7  and  the  analytical 
re.Mill.s  for  i/  — *  oo.  This  diagram  indicates  that,  a 
reduction  in  the  overall  liquid  moments  for  a  gi'~n 
fluid  mass  can  be  achieved  by  splitting  the  cylindri¬ 
cal  volume  into  slices  of  low  aspect  ratio. 

Discussion 

The  codes  previously  in  use  may  serve  for 
establishing  some  basic  results  but  are  too  inefficient 
(and  insufficiently  verified)  for  routine  applications. 
The  finite-difference  code  developed  at  Sandia 
Laboratories  3- 4  rests  on  Chorin’s  method  of  artificial 
compressibility  and  provides  the  steady  solution  at 
11  X  24  X21  grid  points  in  r  ,  <f> ,  z -direction  by 
integrating  over  typically  101  to  8T0*  time  steps,  a 
task  that  requires  G  to  48  minutes  of  CPU  time  on  a 
Cray-lS.  The  result  consists  of  over  22,000  values  of 
the  velocities  v, ,  v  +  ,  v ,  and  the  pressure  p  . 

Strikwerda  &  Nagel  5  briefly  describe  a  code 
using  finite  differences  in  radial  and  axial  direction 
and  pseudospectral  differencing  in  the  azimuthal 
direction.  Nonuniform  grids  are  introduced  for 
increased  resolution  near  the  walls.  The  difference 
equations  are  solved  by  an  iterative  method  based  on 
successive  over-relaxation.  The  computer  time 
required  is  comparable  to  that  of  the  Sandia  code.  A 
thorough  evaluation  of  the  two  codes  is  currently 
conducted  at  BRL.7 

The  experience  with  the  present  version  of  the 
spectral  code  shows  that  high  performance  can  be 
achieved.  The  solution  is  obtained  in  semi-analytica! 
form  with  only  N  ~  4-K  L-M  (typically  less  than 
500)  numerical  coefficients.  This  low  data  volume  is 
especially  attractive  for  storage  and  for  communica¬ 
tion  witli  remote  supercomputers.  The  code  is  very 
well  suited  for  vectorization,  since  practically  all 
CPU  time  is  spent  on  constructing  and  solving  an 
algebraic  system.  The  code  demands  larger  memory 
than  other  codes,  because  (M-hit.  arithmetic  is  highly 
recommended  for  spectral  methods  in  general,  and 
the  algebraic  system  requires  N  (N  +  1)  words  of 
storage.  A  run  with  N  —  500  requires  about  2 
Mbyte  of  memory  and  can  easily  be  carried  out  on 
nowadays  engineering  workstations  within  a  few 
minutes.  (Most  of  our  work  was  done  on  an 
MC68010-based  Apollo  DN300  workstation.  Runs 
with  N  —  500  require  160  minutes,  while  moment 
calculations  with  N  =  192  require  less  than  3 
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minutes  per  poliii.)  Since  the  memory  requirement  is 
acceptable  even  if  higher  resolution  is  desired,  the 
method  applied  here  is  a  viable  alternative  in 
itumcroim  other  lluiil  inn  liimicnl  1 1 r* > l » 1  ••  n i n _  Tin*  abil¬ 
ity  to  obtain  accurate  solutions  for  the  steady  prob¬ 
lem  directly  from  (large)  algebraic  system  bears  valu¬ 
able  potential  to  answer  the  question  whether  the 
steady  solution  is  stable,  and  allows  analysis  of 
unsteady  motions  with  implicit  time-stepping.  The 
design  of  a  reliable  code  for  the  unsteady  problem 
can  take  profit  from  the  knowledge  of  the  eigenvalue 
spectrum  for  small  unsteady  disturbances  of  the 
steady  flow. 

While  the  calculation  of  velocity  and  pressure 
fields  provides  insight  into  the  physics  of  the  flow, 
the  practical  interest  in  the  moments  for  the  quasi- 
steadily  changing  parameters  in  flight  simulations 
can  he  satisfied  with  modest  amounts  of  computer 
time.  This  is  due  to  using  the  modified  Newton 
method  which  updates  the  Jacobian  only  when 
demanded  by  deteriorating  convergence. 

In  general,  the  volume  approach  provides  much 
more  accurate  results  than  the  surface  approach. 
This  is  due  to  the  additional  smoothing  of  fluctuating 
data  by  integrating  over  three  instead  of  two  space 
directions  and  to  using  fewer,  less  fluctuating,  and 
more  accurate  input,  data.  The  absence  of  vr  in  the 
volume  formulation  is  welcome.  This  velocity  com¬ 
ponent  is  small  over  most  of  the  cylinder  length  but 
oscillatory  in  the  radial  direction  (Herbert  et  al. 
1 9S7 )  with  considerable  gradients  near  the  wall.  Near 
the  end  walls,  vr  is  of  the  same  order  as  v ,  with 
steep  gradients  toward  the  end  wall.  Inspection  of 
the  velocity  plots  of  Vaughn  et  el. 4  indicates  that 
these  gradients  were  difficult  to  rtoolve  by  the  finite 
differ-nce  method.  The  aperiodic  component  of  v #  is 
a  p'l, '.lively  small  streaming  term  of  smooth  and 
almost  uniform  i  '-havior  along  the  cylinder  axis. 
The  large  azimuthally  periodic  components  of  vi 
near  the  end  walls  do  not  affect  the  moment  calcula¬ 
tion. 

Probably  the  greatest  advantage  of  the  volume 
formulation  is  the  absence  of  the  pressure  from  the 
moment  equations  This  property  favors  the  use  of 
pressure-free  sets  of  basic  equations,  e.  g.  in  terms  of 
vorticity  or  vector  potential.  The  smaller  number  of 
dependent  variables  can  be  exploited  for  further 
increasing  the  efficiency  liven  in  natural  variable 
formulations,  the  pressure  is  difficult  to  obtain  with 
high  accuracy  because  of  the  invalidity  of  the  equa¬ 
tions  in  the  corners  joining  the  flat  end  walls  to  the 
■•vlindrmal  side  wall  As  shown  in  Figure  3,  the  pres¬ 
sure  may  as.iiirne  extrema  and  vary  strongly  near  the 


corners  and,  therefore,  inaccuracies  in  this  region 
may  influence  yaw  and  pitch  moments.  In  this  con¬ 
text,  it  is  instructive  to  evaluate  the  convergence  his¬ 
tory  uf  t.lir  artificially  tinin-iltipniulmit  lunlliiiil  iltipln 
nieiil, rd  in  the  Snndin  code.3  While  the  velocity 
rapidly  reaches  a  quasi-steady  state,  about  75%  of 

I  lie  il  eral.inns  arc  spent  on  improving  the  pressure 
field'.  We  estimate  that  by  use  of  i.iie  volume 
approach  equivalent  or  superior  values  for  the 
moments  could  be  obtained  with  less  than  20%  of 
the  iterations.  It  is  worthwhile  to  note  that  the 
analytical  results  of  Roscnblat  et  al.,®  and  equations 
(21)  for  the  moments  are  valid  for  closed  containers 
of  more  genera!  shape  and  thus  can  be  used  for  other 
interior  flow  problems. 
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Figure  7.  Contour  plot  of  »he  pressure  field  in  the  Figure  8  Contour  plot  of  the  pressure  filled  in  the 

planes  <t>  =  0'  (left)  and  and  —  90°  (right)  at  plane  s  —  0  25  at  Re  =  300.  Levels  every  0  005. 

Re  ~  300  for  i  >  0.  Levels  every  0.005. 


1  gure  9  Roll  moment  M,  vs  Reynolds  number  Re  for 
,)  —  4.308,  r  —  0.10007.  and  0  -----  20°  Comparison  with 
numerical  and  experimental  data 
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Figure  10.  Yaw  moment  M,  and  pitch 
moment  vs.  Reynolds  number  Re 
for  T)  =  4.368,  r  0.16667,  and 
6  =  20* .  Comparison  with  numerical 
data. 


Figure  11.  Ratios  of  nonlinear  to  linear 
solutions  for  the  moments  at  Re  —  20, 
n  =  4.368. 


Figure  12.  Yaw  moment  M,  and  pitch 
moment  per  unit  length  vs.  aspect 
ratio  r;  at  Re  =  10,  r  =  0.123,  and 
9  =»  2* .  Comparison  of  present  numeri¬ 
cal  results  with  analytical  results  for 
t)  — *  oo  and  data  obtained  by  Nusca 
with  Strikwerda’s  code. 
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ABSTRACT 

The  fluid  motion  in  a  liquid-filled  payload  container  can  jeopardize  the  stable  flight  of  spin- 
stabilized  projectiles.  We  have  developed  analytical  solutions  to  obtain  estimates  and  an 
efficient  numerical  (spectral)  method  to  calculate  accurate  values  of  the  moments  exerted  by 
the  liquid  on  the  projectile.  These  moments  have  been  incorporated  into  a  modified  version  of 
the  numerical  flight  simulator  suggested  by  Vaughn,  Oberkampf  &  Wolfe  Prototype  flights 
with  liquid  versus  solid  payload  veiify  the  destabilizing  effect  of  the  liquid  payload  and  are 
consistent  with  field  tests.  Moreover,  we  demonstrate  that  such  flight  simulations  can  be 
expedited  on  today’s  engineering  work  stations. 

1.  Introduction 

Well-designed  spin-stabilized  projectiles  can  exhibit  a  severe  flight  instability  if  they  carry 
a  liquid  payload.  Two  types  of  such  instability  can  be  excited  by  the  coning  motion  of  the 
projectile  about  the  trajectory.1  The  first  type  originates  from  resonance  with  inertial  waves 
at  critical  frequencies  (ratios  of  coning  rate  Q  to  spin  rate  w).  This  resonance  is  most  pro¬ 
nounced  for  low-viscosity  liquid  fills,  i.e  at  high  Reynolds  numbers,  and  depends  sensitively  on 
the  cylinder’s  aspect  ratio.  Theoretical  analysis  involving  the  boundary-layer  approximation 
provides  design  criteria  for  sufficiently  large  Reynolds  numbers,  It  >  1000,  say.  We  define  the 
Reynolds  number  by  R  —  pu>a2/p,  where  p  is  the  density,  a  the  radius  of  the  cylindrical 
cavity,  and  p  the  viscosity.  The  second  type  of  instability  is  of  essentially  viscous  nature  and 
occurs  at  iow  and  medium  Reynolds  numbers  for  a  wide  range  of  aspect  ratios  and  frequencies. 
Theoretical  analysis  must  be  based  on  the  Navier-Stokes  equations.  For  a  wide  range  of  Rey¬ 
nolds  numbers,  both  types  of  instability  may  appear  simultaneously. 

To  enable  the  design  of  reliable  projectiles,  the  effects  of  the  liquid  payload  must  be 
analyzed  and  incorporated  into  the  design  tools  such  as  the  numerical  flight  simulators.  A 
prototype  flight  simulator  (FFSfiDOF)  has  been  developed  by  Vaughn  ct  al.2  by  integrating  a 
table  of  liquid  moments  into  a  six-degiee-of-freedom  code.  The  moments  were  computed  with 
a  Navier-Stokes  solver  (FFSG)3  '1  The  computational  demands  of  this  two-level  code  prevent 
application  in  practice. 

In  the  following  we  describe  two  more  efficient  procedures  to  obtain  the  liquid  moments 
for  integration  into  flight  simulations  The  first  method  rests  on  quick  estimates  of  the 
moments  from  analytical  solutions,  the  second  efficiently  calculates  accurate  moments  in  a 
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wide  range  of  Reynolds  numbers  by  means  of  a  spectral  Navier-Stokes  code.  The  development 
of  a  new  method  for  calculating  the  moments  was  essential  in  achieving  these  goals.  One  of 
the  flights  studied  bv  V'aughn  et.  al  is  used  as  an  example  to  show  some  shortcomings  of  the 
earlier  approach  and  to  demonstrate  the  feasibility  of  flight  simulations  for  liquid-filled  shells 
on  engineering  work  stations. 

2.  Approximate  Solutions 

In  earlier  work,5  6  we  have  reported  results  of  a  perturbation  analysis  of  the  nonlinear 
problem  for  the  deviation  of  the  velocity  from  rigid-body  rotation  in  a  finite  segment  of  an 
infinite  cylinder.  The  expansion  parameter  is  e  =  (ft  /<*.’)  sintf,  where  0  is  the  coning  angle. 
At  the  linear  level,  0(c),  we  obtain  a  purely  axial  flow  in  the  form 
i,i  =  -2c  Im{|/|(or ) Tjfo)  -  r  where  /]  denot  s  the  modified  Bessel  function. 

q  —  (1  -  i  )[R  2)1  and  the  azimuthal  angle  <?  is  measured  from  the  plane  of  the  two  axes 
of  rotation  This  component  is  the  dominating  feature  of  the  flow  in  a  long  cylinder  The 
second  order  terms,  0{(2).  were  also  obtained  in  analytical  form.  Of  these  terms,  only  the 
aperiodic  component  of  the  azimuthal  velocity  is  relevar  t  to  the  moments.  This  component 
turns  out  to  be  proportional  to  (lie  component  of  v.  j  in  the  plane  ©  -  .2. 
v  <?o  =  ~  2c  lle{/ i(ar  )/ j(c»)  -  r  }  .  This  solution  is  valid  for  arbitrary  Reynolds  number. 
Owing  to  the  neglect  of  the  end  walls,  however,  this  analytical  result  is  expected  to  he  a  rea¬ 
sonable  approximation  only  for  cylinders  of  sufficiently  large  aspect  ratio. 

3.  Spectral  Navior  Stokes  Code 

To  calculate  the  deviation  of  velocity  field,  pressure  field,  and  moments  from  solid-body 
motion  in  a  wide  range  of  parameters  for  finite-length  cylinders,  we  have  developed  a  fully 
spectral  code  1  for  solving  the  Navier-Stokes  equations.  The  code  uses  Fourier  series  in  the 
azimuthal  direction  and  combinations  of  Chebyshev  polynomials  in  the  radial  and  axial  direc¬ 
tions  The  combinations  arc  chosen  such  that  the  no-penetration  and  no-slip  boundary  condi¬ 
tions  are  implicitly  satisfied.  By  appropriate  choice  of  even  or  odd  polynomials,  we  avoid 
singularities  on  the  axis  and  implement  the  symmetries  of  the  solution  in  the  upper  and  lower 
half  cylinder.  Gauss-Lobatto  points  are  used  for  collocation  such  that  no  points  are  locat'd  or: 
.it1  surface  nor  the  axis.  This  choice  avoids  the  occurrence  of  spurious  pressure  terms  associ¬ 
ated  with  the  corners  of  the  domain  and  solves  the  axis  problem  in  an  elegant  way. 

The  solution  procedure  consists  of  two  basic  elements,  the  first  of  which  rests  on  a  linear¬ 
ization  of  the  Navier-Stokes  equations  about  some  known  solution  vrf  .  This  ‘Known’  solution 
may  be  =  0,  or  in  general,  a  solution  previously  obtained  for  neighboring  parameters 
Tiie  collocation  method  converts  the  linear  partial  differential  equations  into  a  linear  algebraic 
system  of  dimension  ,Y  —  4  K  L  M  for  the  coefficients  of  the  natural  variables  v/,  v %  ,  r/. 
and  p  4 .  The  numbers  A',  L,  At  of  expansion  functions  in  radial,  azimuthal,  and  axial  direc¬ 
tion,  respectively,  govern  the  accuracy  of  the  approximation. 

The  second  element,  of  the  procedure  utilizes  the  modified  Newton  method  for  solving  the 
nonlinear  problem.  In  the  modified  form,  the  Jacobian  is  not  updated  in  every  step  but  only 
when  necessary.  The  need  to  recalculate  th**  J?cobi»n  is  detected  by  a  simple  algorithm  that 
monitors  the  rate  of  convergence.  This  strategy  is  especially  suited  lor  flight  simulations 
where  the  parameters  vary  slowly  with  the  flight  time. 
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4.  Volume  Integrals  for  the  Moments 

In  earlier  work,3'  '*•  8'  9  the  moments  are  calculated  by  integration  of  normal  and  tangen¬ 
tial  stresses  acting  at  the  inside  wall  of  the  cylinder  This  ‘surface  approach’  derives  the 
moments  from  the  pressure  and  the  velocity  gradients  at  the  wall.  These  values  are  difficult  to 
obtain  with  good  accuracy.  Li  &  Herbert  10  have  developed  a  new  formulation  for  the 
moments  which  rests  on  volume  integrals.  The  nondimensional  moment  coefficients  can  be 
shown  to  take  the  form 

ij  2  it  1 

=  — — — J  J  J  v,  r  2eoso  drd  $  dz  ,  A//—-  A// tan  9  , 
tan®  .  ,  o  o 

i|  !i  I  i)  2*  1 

=  ef  J  f  v ^r~ drd  <t> dz  -i - j  f  J  t'*r2sin d>drd<f>dz  , 

-  i)  o  0  t,an"  -  i)  o  o 

where  r/  is  the  cylinder’s  aspect  ratio.  The  reference  moment  is  p^2as.  The  volume  integral 
approach  involves  only  the  radial  and  azimuthal  velocity  components  Integration  over  do 
reduces  the  requirements  to  the  knowledge  of  the  aperiodic  component  t’^0  of  e 9  and  the  sim¬ 
ply  periodic  component  t>,j  of  v2.  Therefore,  the  volume  approacn  can  also  be  applied  to  the 
analytical  results  given  above  and  provides  yaw  and  pitch  moments  without  explicit  knowledge 
of  the  pressure 

5.  Results  for  the  Liquid  Moments 

While  the  accurate  descr  IpkiCjj  uf  velocity  and  pressure  at  higher  Reynolds  numbers 
requires  many  expansion  functions,  the  smoothing  of  small  oscillations  by  the  volume  approach 
allows  calculation  of  the  moments  from  crude  approximations  The  (absolute)  accuracy  of 
10' 3  for  engineering  applications  can  be  achieved  with  spectral  approximations  as  low  as 
I\  =  q,  L  =  3,  A/  =  4.  The  accuracy  has  to  be  seen  in  the  light  of  considerable  uncertainty 
in  the  moments  governing  the  exterior  aerodynamics  of  the  projectile.  We  note,  however,  that 
large  aspect  ratios  may  require  additional  expansion  functions  in  axial  direction  while  increas¬ 
ing  Reynolds  number  in  general  requires  higher  resolution  in  both  radial  and  axial  direction. 

Figure  1  shows  the  comparison  of  numerical  results  for  the  yaw  and  pitch  moments.  The 
roll  moment  is  proportional  to  the  yaw  moment.  While  the  agreement  for  the  yaw  moment  at 
high  Reynolds  numbers  is  surprisingly  good,  the  deviation  in  the  pitch  moment  is  likely  to  ori¬ 
ginate  from  insufficient  resolution  of  the  steep  pressure  gradients  in  the  finite-element  code  9 
The  effect  of  discretization  errors  has  been  reduced  in  the  results  of  the  spectral  code  by  use  of 
the  volume  integrals  for  the  moments.  Figure  2  compares  the  analytical  estimates  for  yaw  and 
pitch  moments  with  the  numerical  results  of  the  spectral  code.  Some  systematic  deviations 
exist,  hut.  t.he  quantitative  and  qualitative  features  are  very  similar  The  computational  effort 
for  the  analytical  results  is  negligible. 

0.  A  Prototype  Flight 

We  have  used  the  flight  of  an  M483  projectile  as  described  by  Vaughn  et  a!  2  as  a  test 
case  for  verifying  the  function  of  our  version  of  the  flight  simulator,  for  studying  the  effect  of 
th**  improved  moments  in  comparison  with  earlier  work,  and  for  evaluating  the  different  con¬ 
cepts  of  incorporating  the  fluid  moments.  Some  erratic  behavior  of  the  original  code 
FFS6DOF  under  extreme  flight  conditions  will  be  cultivated  in  a  future  version. 
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Two  options  have  been  considered  for  incorporating  the  liquid  moments.  The  first  way 
utilizes  a  table  of  moments  and  linear  interpolation,  as  the  original  code.  Alternatively,  the 
moments  are  calculated  by  a  subroutine  which  is  called  whenever  the  parameters  vary  beyond 
some  threshold.  The  latter  approach  avoids  the  error  introduced  by  linear  interpolation  on  a 
coarse  grid  and  is  especially  efficient  for  stable  flights  after  an  initial  transient  period  after 
launch.  For  systematic  studies  of  different  launch  conditions  for  a  given  shell,  however,  the 
precalculated  table  seems  more  appropriate  and  the  efficiency  of  the  spectral  code  permits  use 
of  finer  grids  Some  accuracy  can  aLo  be  gained  by  improving  the  interpolation. 

Our  computational  work  has  been  performed  on  a  Sun  3/1-10  work  station  with  floating 
point  accelerator.  Simulation  of  a  typical  flight  (32  seconds  flight  time)  requires  12  minutes  if 
the  liquid-moment  table  is  precalculated  Generation  of  this  table  from  the  analytical  solution 
requires  less  than  10  seconds.  Incorporation  of  the  analytical  estimates  by  subroutine  contri¬ 
butes  a  negligible  increase  in  the  simulation  time.  Precalculation  of  the  9x8x8  moment 
table  with  K  —  4,  L  —  3,  and  A 1  =  4  requires  less  than  2  hours,  although  the  sequence  of 
the  calculations  has  not  yet  been  optimized.  This  sequence  is  relevant  when  using  the 
modified  Newton  method 

In  the  following  cases,  we  keep  all  parameters  at  the  values  specified  by  Vaughn  et  al  2 
except  the  pitch  rate  at  launch  The  initial  data  for  the  flight  simulation  -  angle  and  pitch  rate 
:»t.  tt-p  mu^-df*  -  v-e  d-ffi<’n]r.  fo  determine  and  may  vary  between  flights.  An  initial  pitch  rate 
of  -1.9-16  rad/s  was  chosen  by  Vaughn  et  al  to  match  the  observed  flight  behavior, ^  in  partic¬ 
ular  the  truncated  flight  time  of  approximately  26  seconds  Figure  3a  shows  the  history  of  the 
aeroballistir  angle  of  attack  for  a  solid  payload  (with  all  liquid  moments  set  to  zero)  that 
represents  the  ‘normal’  flight  of  a  projectile  of  this  type.  The  results  of  Vaughn  et  al  (ref.  2, 
figure  3)  for  a  liquid  payload  are  reproduced  with  our  version  of  the  code  in  Figure  3b.  Figure 
3c  shows  the  result  for  the  same  initial  conditions  with  the  liquid  moments  calculated  by  the 
spectral  code.  Use  of  the  analytical  estimates  leads  to  an  almost  identical  picture.  With  the 
correct  moments,  the  flight  is  essentially  stable,  with  some  high-frequency  motion  over  the 
whole  flight  time.  The  liquid  payload  does  not  prevent  the  shell  from  achieving  the  full  flight 
time  or  distance. 

When  using  the  correct  moments,  we  observe  the  flight  time  reported  by  Vaughn  et  al.2 
and  observed  in  the  field  tests  11  with  initial  pilch  rates  in  excess  of  -3.5  rad/s  as  shown  in 
Figure  4a  Under  the  given  transonic  launch  conditions,  this  pitch  rate  appears  very  realistic 
(D'Amico,  personal  communication).  Simulation  of  the  same  flight  using  the  analytical  esti¬ 
mates  for  the  moments  leads  to  earlier  onset  of  the  instability,  as  shown  in  Figure  4b.  A  flight 
time  of  26  seconds  is  obtained  with  the  lower  initial  pitch  rate  of -3  2  rad  s.  Since  the  usp  of 
the  analytical  estimates  predicts  less  stable  behavior,  the  estimates  can  be  used  for  a  first  con¬ 
servative  check  of  the  design. 


7.  Summary 

We  have  developed  analytical  solutions  lor  the  liquid  moments  which  aiiow  quick  and 
conservative  estimates  on  the  flight  instability  of  liquid-filled  shells.  For  more  accurate  stu¬ 
dies,  the  liquid  moments  can  be  generated  with  a  spectral  Navier-Stokes  solver  for  Reynolds 
numbers  in  the  range  R  <  1000  While  the  analytical  results  predict  only  the  viscous  flight 
instability,  the  Navier-Stokes  solutions  incorporate  b:th, 

1 1- v  due  to  resonance  with  inertial  waves 


H?,rJ 


1  instabil- 
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The  computation  of  the  table  of  liquid  moments  in  the  most  critical  range  of  Reynolds 
numbers,  R  <  100,  requires  less  than  two  hours  on  an  engineering  work  station  Sun  3/110 
FPA.  The  flight  simulation  runs  at  approximately  25  times  real  time,  typically  12  minutes 
These  data  clearly  show  that  the  liquid  moments  can  be  incorporated  into  the  flight  simula¬ 
tions  for  practical  applications. 
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Figure  3  Aeroballistic  angle  of  attack  as  a  function  of  the  time  of  flight  Initial  pitch  rate 
-1  9-lC  rad/s.  (a)  Solid  payload,  (b)  Liquid  payload.  Table  of  liquid  moments  produced  by  the 
code  FFSG  of  Vaughn  et  al  (c)  Liquid  payload.  Tabic  of  liquid  moments  from  spectral  code. 
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ABSTRACT 


Artillery  shells  with  liquid  payloads  may  experience  a  severe  flight  instability 
owing  to  the  moments  exerted  by  the  viscous  fluid  motion  in  the  cylindrical  pay- 
load  container.  Incorporation  of  these  moments  into  flight  simulations  as  a  rou¬ 
tine  design  tool  requires  a  highly  efficient  code  for  solving  the  Navier-Stokes 
equations.  We  describe  a  spectra!  collocation  method  which  is  based  on 
Chebyshev-Fourier-Chebyshev  expansions  in  the  radial,  azimuthal,  and  axial 
direction.  The  method  exploits  the  symmetries  of  the  problem.  Using  a  volume 
approach  and  an  analytical  result  by  Rosenblat,  accurate  moments  are  obtained 
in  small  fractions  of  the  time  required  by  other  codes.  Solutions  for  the  steady 
motion  are  presented  and  compared  with  numerical  and  experimental  data. 

Introduction 

Gyros  and  rotating  fluids  often  exhibit  unexpected  behavior.  In  the  past,  it 
has  been  recognized  that  spin-stabilized  shells  with  liquid  payloads  can  suffer  a 
dynamical  instability  originating  from  resonance  with  inertial  waves.1  Since  this 
phenomenon  is  basically  inviscid  and  is  routinely  avoided  by  proper  design,  it  was 
surprising  to  observe  in  some  cases  another  type  of  instability  which  is  character¬ 
ized  by  an  increase  in  nutation  (or  yaw)  angle  and  a  simultaneous  loss  in  spin 
rate.  The  rapid  drop  in  spin  rate  is  clearly  a  viscous  phenomenon,  and  labora¬ 
tory  experiments,  computational  results,  and  field  tests  have  meanwhile  shown 
that  this  instability  is  caused  by  the  nutation-induced  fluid  motion  in  a  certain 
range  of  relatively  small  Reynolds  numbers.  Although  in  special  cases  this  insta¬ 
bility  has  been  overcome  by  trial  and  error,  future  design  of  reliable  projectiles 
would  take  profit  from  the  opportunity  to  calculate  the  liquid  moments  and  to 
account  for  these  moments  in  flight  simulations.  The  empirical  data  base  is 
sparse,  however,  and  the  computational  methods  in  use  3> 6  are  rather 
demanding.  An  evaluation  and  verification  of  the  codes  by  1 'aughn  et  ah'4  and 
Strikwcrda  &  Nagel  ,c>  is  currently  conducted  at  BRL.^  Typical  computer  times 
for  a  single  case  are  in  the  range  of  G-12  hours  on  VAX-class  machines.  Six- 
degree-of-freedorn  flight  simulations  8  typically  use  2  105  time  steps  over  the 
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flight  time  of  the  order  of  30  seconds.  Study  of  the  interaction  of  the  interior 
fluid  motion  with  the  exterior  aeroballistics  consequently  requires  either  a  very 
fast  subroutine  for  calculating  the  liquid  moments  or  interpolation  in  a  multi¬ 
dimensional  table  of  500  -  1000  8  precalculated  values.  Hence,  flight  simulations 
for  liquid-filled  shells  are  currently  a  very  expensive  tool  and  are  not  ready  for 
routine  applications. 

In  previous  work,9  we  conducted  a  theoretical  analysis  which  aimed  at  the 
origin  of  the  viscous  despin  (negative  roll)  moment  in  cylinders  of  large  aspect 
ratio.  This  analysis  showed  that  the  deviation  from  solid-body  rotation  is 
governed  by  a  small  parameter,  c  =  (0  /w)sin0,  involving  the  nutation  rate  fl  , 
the  nutation  angle  6,  and  the  spin  rate  u>.  A  solution  of  the  linearized  equations 
was  developed  for  a  finite-length  segment  of  an  infinitely  long  cylinder,  i.  e.  disre¬ 
garding  the  end  walls  of  the  cylinder.  Velocity  field  and  the  viscous  components 
of  the  moments  were  obtained  in  closed  form.  The  velocity  field  consists  only  of 
an  axial  component  of  order  O  (e)  which  is  the  prominent  feature  of  the  fluid 
motion  in  slender  cylinders  and  produces  a  negative  roll  moment  of  order  0(e2) 
owing  to  Coriolis  forces.  Although  this  roll  moment  is  in  reasonable  agreement 
with  experimental  and  computational  data,  the  original  analysis  accounted  only 
for  the  viscous  part  of  the  yaw  and  pitch  moments.  These  latter  moments  con¬ 
tain  essential  contributions  of  the  pressure  4  that  originate  from  the  turning  of 
the  flow  near  the  end  walls  and  were  not  captured  by  the  linear  analysis.  The 
effect  of  nonlinearity  was  studied  10  by  using  perturbation  expansions  in  t  and 
was  found  to  be  small  except  for  an  aperiodic  streaming  term  in  the  azimuthal 
direction. 

Whereas  the  pressure  field  cannot  be  captured  by  the  perturbation  approach, 
it  provided  valuable  insight  into  the  structure  of  equations  and  solution.  The 
analytical  work  suggests  the  use  of  a  numerical  method  that  exploits  (i)  the 
near-linearity  of  the  governing  equations  and  (ii)  the  smoothness  of  the  solution 
in  the  relevant  range  of  Reynolds  numbers.  We  have  therefore  pursued  a  simple 
concept  that  is  open  to  further  refinements.  We  use  Chebyshev-Fourier- 
Chebyshev  expansions  in  respectively,  and  convert  the  linearized 
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equations  into  an  algebraic  system  for  the  expansion  coefficients.  Linearization 
can  be  performed  about  the  trivial  solution  or  any  other  known  solution,  e.  g.  at 
neighboring  parameters.  The  solution  of  the  linear  algebraic  system  is  used  as 
initial  approximation  for  iterative  improvement  by  the  modified  Newton  method. 
The  feasibility  of  this  approach  has  been  demonstrated  11  with  a  crude  spectral 
approximation  to  the  solution.  Problems  in  calculating  the  pressure  that  arise 
from  the  invalidity  of  the  basic  equations  along  the  joint  of  the  flat  end  walls  to 
the  cylindrical  side  wall  have  meanwhile  beer,  overcome.12  The  present  version  of 
the  code  exploits  the  diametral  symmetry’  of  the  flow  about  the  center  of  the 
cylinder  and  allows  for  higher  resolution  at  modest  CPU  times.  This  version  can 
also  be  adapted  lor  the  analysis  of  unsteady  problems.  Dramatic  increase  in 
efficiency  has  recently  been  achieved  13  by  combining  an  analytical  result  of 
Rosenblat  et  al.6  with  a  volume  formication  for  calculating  the  liquid  moments. 
The  moments  can  be  obtained  from  only  the  simply  periodic  components  of  the 
axial  velocity  and  the  azimuthal  streaming  term.  A  fast  subroutine  for  flight 
simulations  exploits  the  analytical  results.  For  more  accurate  studies,  complete 
tables  of  moments  can  be  calculated  in  a  few  hours  on  an  engineering  worksta¬ 
tion. 


Governing  Equations 

We  consider  the  steady  motion  of  a  fluid  of  density  p  and  viscosity  p  in  a 
cylinder  of  radius  a  and  length  2c  in  an  aeroballistic  coordinate  system  x  ,y ,  z , 
where  z  is  the  axis  of  the  cylinder,  as  shown  in  Figure  1.  The  inertia!  axis  Z  in 
flight  direction  and  the  2 -axis  enclose  the  nutation  angle  9.  The  cylinder  rotates 
with  the  spin  rate  u  about  z  while  the  x,  2 -plane  rotates  with  the  nutation  rate 
0  about  the  if -axis.  Spin  rate  u>  and  nutation  rate  0  are  constant.  All  quanti¬ 
ties  are  made  nondimensional  using  a  ,  w,  and  p  for  scaling  length,  time,  and 
mass,  respectively.  The  solution  depends  on  four  nondin  ensional  parameters: 
aspect  ratio  75  =  c  /a  ,  nutation  angle  9,  frequency  r  =  0  /w,  and  Reynolds 
number  Re  —  puja2fp.  The  aspect  ratio  enters  the  solution  only  through  the 
boundary  conditions  at  the  end  walls  of  the  cylinder.  The  motion  is  subject  to 
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the  nosiip  and  no-penetration  conditions  at  the  cylinder  walls.  Since  the  velocity 
field  degenerates  for  either  =  0,  Cl  —  0,  0  =  0,  or  p—*oo  to  rigid-body  rota¬ 
tion  of  the  fluid,  it  is  appropriate  to  concentrate  on  the  deviation  vd  of  the  velo¬ 
city  from  rigid-body  rotation, 

V  =  vr  +  \d  ,  vf  =  r  e0  (la,b) 

where  is  the  azimuthal  unit  vector.  The  boundary  conditions  on  vd  are 
homogeneous.  The  pressure  field  is  split  according  to 

p  =  pr  +  pd  ,  pr  =  r2(l+rj2  +  r%2  +  *2e2  -  2  rz  n  T  (2a, b) 


where  rf  =  -ecos <f>,  t ^  =  esin<£,  tz  =  rcos0,  e  =  rsin#.  The  pressure  pr 
differs  from  the  pressure  in  rigid-body  rotation.  The  form  of  p r  is  chosen  such 
that  the  reduced  pressure  pd  appears  only  in  the  z  -momentum  equation. 

In  cylindrical  coordinates  r  ,<t>,  z ,  the  equations  for  the  velocity  components 
\d  —  (t’r ,  Vj,  vz )  and  pressure  p  d  take  the  form 


1  a ,  v  ,  i  dv*  ,  dvt 

— — (rur)+ - —  +  -t—  =  0 

r  or  r  o4>  dz 


(3a) 
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The  primary  effect  of  nutation  is  the  ^-periodic  force  term  -2 rrf  =  2ercos<£  in 
the  z  -momentum  equation  (3d).  For  e  =  0,  equations  (3)  have  the  trivial  solu¬ 
tion  \i  ~  0,  p  d  =  0.  The  system  supports  the  following  symmetries: 

vr(r,  <p+n,-z)  =  vr(r  ,  <j> ,  z)  ,  v  ^(r  ,  <t>+n ,  -  z  )  =  t/0(r  ,  <f>,  z  )  (4a, b) 

v:(r  ,  4>+tt  ,  -  z)  =  -  v2(r  ,  4>,  z)  ,  pd  (r  ,  <j>+n  ,  -  z)  =  pd  (r  ,  <j>,  z)  (4c, d) 

Therefore  it  is  sufficient  to  obtain  a  solution  in  the  half-cylinder,  z  >  0,  with 
appropriate  symmetry  conditions  at  z  =  0. 


Some  Analytical  Results 

The  steady  flow  in  a  relatively  long  cylinder  (aspect  ratio  rj  >  4)  at  low 
Reynolds  number  is  expected  to  exhibit  little  axial  variation  over  most  of  the 
cylinder  length.  Previous  work  9  has  therefore  relaxed  the  boundary  conditions 
at  the  end  walls  and  studied  the  steady  flow  in  a  finite  segment  of  an  infinitely 
long  cylinder. 

In  the  physical  situations  of  interest,  e  =  (Cl  /u)  sin#  is  a  small  parameter, 
and  consequently,  it  is  reasonable  to  pursue  a  perturbation  expansion  in  e.  This 
provides  \d  in  the  form 

v'  =  £e"v<"  >(,.,*)  (5) 

n  =  1 

and  similar  expressions  for  pd ,  The  development  of  expressions  for  the  expan¬ 
sion  coefficients  v(n'  from  equations  (3)  leads  to  an  alternating  pattern: 


(0,  0,  ,  n  odd 

(vr(n),  Vj  (n\  0)  ,  n  even 


t<5) 

(7  a) 
(7b) 
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(7c) 


„,(» )  J"  e  «>»  -  ■)♦) 

m  =  1 


where  the  tilde  denotes  the  complex  conjugate.  The  aperiodic  term  in  vf^”^  is 
suppressed  by  the  continuity  equation.  The  r  -dependent  coefficient  functions  in 
Eqs.  (7)  are  required  to  satisfy  homogeneous  boundary  conditions  at  r  —  1  and 
to  be  finite  at  the  axis  r  =  0. 

The  axial  velocity  at  order  0(f)  can  be  found  in  analytical  form, 

7,(ar  ) 

"1.(0=  •  l  ,  ,  '  -  0  (8) 

/l(a) 

where  I  x  denotes  the  modified  Bessel  function,  and  a  =  (1  +  i)(Rc  /2)1/2.  This 
solution  is  valid  for  arbitrary  Reynolds  number  but  may  be  unstable  as  Re 
exceeds  some  critical  value.  This  component  is  the  dominating  feature  of  the 
flow  in  a  long  cylinder.  The  interesting  properties  of  the  associated  flow  field  are 
discussed  by  Herbert.8 

At  order  O  (e2),  comparison  of  the  equation  for  v2o  with  the  imaginary  part 
of  the  equation  for  u/n  immediately  shows  that  the  aperiodic  component  of  the 
azimuthal  velocity  is 

”20  (0=  —  2Im[u/n(r)1  (9) 

The  ^-periodic  components  are  governed  by  a  coupled  set  of  inhomogeneous 
differential  equations  with  variable  coefficients.  With  some  effort,  the  radial  velo¬ 
city  component  of  0(e2)  can  be  found  in  closed  form.10  In  view  of  the  effort 
involved  in  deriving  the  analytical  result  and  the  ultimate  need  to  determine  cer¬ 
tain  coefficients  numerically,  the  differential  equations  for  the  third-order  com¬ 
ponents  were  solved  by  means  of  a  spectral  collocation  method. 

The  motion  is  governed  by  the  axial  component  wn  at  order  0(f)-  Of  the 
higher  order  terms,  only  the  aperiodic  term  v2o  is  substantial.  In  the  cylinder’s 
center  section,  these  terms  are  in  good  agreement  with  computational  results.  All 
the  other  terms  are  not  only  of  order  0(1)  but  in  fact  less  than  unity,  assuring 
rapid  convergence  of  the  perturbation  series.  The  contribution  of  ic3)  to  the 
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despin  moment  is  negligible.  The  (^-periodic  terms  oscillate  about  zero  as  r 
varies  between  0  <  r  <  1.  Accurate  representation  of  single  high-order  terms 
by  radial  Chebyshev  series  may  require  numerous  expansion  functions.  For  the 
total  velocity  field,  however,  the  error  in  representing  these  terms  is  of  little 
importance.  At  Reynolds  numbers  in  the  range  of  maximum  despin  moment, 
reasonably  accurate  approximations  can  be  obtained  with  only  a  few  polynomials 
in  radial  direction.  In  the  azimuthal  direction,  the  solution  is  governed  by  terms 
periodic  in  <f> ,  and  by  the  aperiodic  term  u20-  Fourier  series  with  three  or  five 
modes,  therefore,  provide  approximations  of  sufficient  accuracy  for  practical  pur¬ 
pose.  The  perturbation  analysis  clearly  shows  that  the  main  features  of  the  flow 
are  governed  by  the  linear  0(e)  part  of  equations  (3)  with  small  corrections  for 
nonlinearity.  This  property  will  not  change  for  a  finite-length  cylinder. 

Spectral  Approximations 

The  results  of  the  analytical  work  suggest  that  a  good  approximation  to  the 
flow  in  a  finite  cylinder  can  be  obtained  by  solving  linearized  versions  of  equa¬ 
tions  (3).  Linearization  can  be  performed  in  different  ways.  The  first  is  a  lineari¬ 
zation  in  e.  Besides  Eqs.  (4),  the  resulting  equations  support  the  additional  sym¬ 
metries 

vrf  (r  ,  0+tt,  a  )  =  -  \d(r  ,  <t>,  z)  ,  p  d  (r  ,  <j>+TT  ,  z  )  =  -  p  d  (r  ,  <t> ,  z  )  (lOa,b) 

These  relations  permit  useful  checks  on  the  results  of  the  spectral  code.  A  second 
linear  system  can  be  obtained  by  linearization  in  the  components  of  \d  .  This 
linearization  retains  coupling  terms  such  as  1ri)vz  in  Eq.  (3b)  which  destroy  the 
symmetries  (10).  The  second  system  can  be  considered  a  special  case  with 
=  0  of  a  linearization  about  some  known  solution  \d .  The  latter  procedure 
is  very  efficient  if  the  solution  is  sought  for  a  densely  spaced  sequence  of  parame¬ 
ter  combinations  as  in  flight  simulations. 

The  algebraic  form  of  the  equations  is  obtained  by  use  of  spectral  colloca¬ 
tion.  The  velocity  components  are  expressed  in  the  form 

<V  =  E  £  E  %.,«,;(’•)  M«!  ^1.(7)  (11) 

*  =  1  1  =  1  rn  =  1  ' 
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with  similar  expressions  for  v vx  ,  and  p  “  .  The  azimuthal  functions  are 

l  -  1 
cos  — 

=  1  J 


• <j>  ,  /  odd 


sin  — d>  ,  l  even 
2 


The  azimuthal  collocation  points  are  equidistant, 

<h  =  2tt  (/  -  1)/L  ,  f  =  1,2, 


(12) 


(13) 


and  L  is  odd. 

In  a  first  version  of  the  code,  radial  and  axial  collocation  points  are  located 
at  the  maxima  of  the  highest  Chebyshev  polynomials.  The  boundary  conditions 
are  implemented  by  replacing  three  of  the  four  differential  equations  in  the  boun¬ 
dary  points.  The  question  then  is  which  equation  should  be  retained  and  where 
the  condition  on  the  pressure,  e.g.  p  d  =  0,  should  be  applied.  Trial-and-error 
leads  to  numerous  cases  with  ill-determined  matrices  or  zero  determinant.  In 
other  cases,  a  correct  solution  for  the  velocity  field  is  obtained,  but  the  pressure 
contains  a  non-physical  spurious  term.  Problems  with  spectral  calculations  of  the 
pressure  in  closed  domains  with  corners  are  well-known  but  the  reports  on  their 
origin  and  methods  for  solution  are  rather  unspecific.  We  have  therefore  per¬ 
formed  a  detailed  analysis  of  the  flow  in  a  square  driven  by  an  internal  force 
field.  This  simpler  two-dimensional  problem  exhibits  all  characteristics  -  includ¬ 
ing  the  spurious  pressure  term  -  of  the  original  problem.  The  study  reveals  that 
the  spurious  term  vanishes  in  all  collocation  points  except  the  corners,  where  it 
may  assume  arbitrary  values.  The  term  can  be  suppressed  by  retaining  in  the 
corners  one  of  the  momentum  equations  that  contain  the  derivative  of  the  pres¬ 
sure  in  the  direction  of  the  boundary. 

In  a  second  version  of  the  spectral  code,  the  problems  of  the  pressure  calcu¬ 
lation  have  been  avoided  by  using  a  different  set  of  collocation  points.  The 
expansion  functions  in  radial  and  axial  direction  depend  on  the  index  /  and  may 
be  different  for  the  variables  vr ,  v#,  vx ,  and  p  d  .  They  are  combinations  of  even 
or  odd  Chebyshev  polynomials  such  that 
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(i)  the  homogeneous  boundary  conditions  are  implicitly  satisfied, 

(ii)  the  symmetry  conditions  (4)  are  satisfied,  and 

(iii)  the  limit  value  of  the  variables  for  r  — ►  0  (i.  e.  the  value  on  the  axis)  is 
independent  of  <j>. 

The  collocation  points  are 

Tk  —  sin  ^  n  ^  =  1,2,  K  (Ida) 

_ZL  =  sin  £ m  ■—  3,  t r  ,  m  =  1,  2,  •  •  •  M  (14b) 

r)  AM 

Since  0  <  rk  ,  no  points  are  located  on  the  axis.  Also,  rk  <  1 ,  zm  <  77 ,  such 

that  no  points  are  located  on  the  surface.  The  points  in  radial  and  axial  direction 

are  concentrated  near  the  boundary  such  that  high  resolution  in  this  region  is 
obtained  without  additional  coordinate  stretching.  Thus  the  boundary  layers 
forming  at  higher  Reynolds  number  can  be  resolved  by  slightly  increasing  K  and 
M. 

The  spectral  collocation  method  converts  the  linear  system  of  partial 
differential  equations  derived  from  Eqs.  (3)  into  an  algebraic  system  of  dimension 
N  —  A  K  L  -M  for  the  coefficients  uktm  ,  vklm  ,  wklm  ,  and  pktm  of  vr ,  v^,  vt ,  and 
p  d  ,  respectively.  The  linear  system  for  the  expansion  coefficients  is  solved  by 
Gauss  elimination  with  partial  pivoting.  The  subroutine  used  retains  all  data 
required  to  solve  the  same  system  with  a  new  right-hand  side  without  repeating 
the  costly  reduction  of  the  matrix  to  upper  triangular  form.  Once  the  solution  is 
obtained,  a  new  right-hand  side  is  formed  taking  the  nonlinear  terms  into 
account  and  the  system  is  iteratively  solved  until  sufficient  accuracy  is  achieved. 
The  procedure  is  equivalent  to  the  modified  Newton  iteration  (without  updating 
the  Jacobian  in  every  step)  and  converges  rapidly  since  the  nonlinear  corrections 
to  the  velocity  are  small  while  the  pressure  appears  linear  in  equations  (3). 
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Results  for  Velocity  and  Pressure 

In  the  following,  we  present  some  results  for  the  velocity  and  pressure  fields 
at  0  =  20°  ,  r  =  0.16667,  and  r\  =  4.368  which  results  in  e  =  0.057.  The 
results  are  for  K  —  6,  L  =  5,  and  M  =■  8,  and  consequently  N  =  960.  Calcu¬ 
lation  of  a  single  solution  with  this  high  resolution  requires  about  2  minutes  on  a 
Cray-lS.  Figure  2  shows  the  axial  and  radial  velocities  in  the  planes  <f>  =  45° 
and  <p  =  135°  at  Re  =  20.  Only  the  upper  half,  z  >  0,  of  the  cylinder  is 
shown;  the  lower  halt  is  governed  by  the  symmetries  (4).  The  scale  values  give 
the  velocity  per  unit  length  where  the  diameter  is  four  units.  The  velocity  distri¬ 
bution  at  z  —  0  agrees  well  with  the  results  of  the  perturbation  analysis  and 
computations  with  the  Sandia  code.3  Near  the  end  walls,  the  solution  is  more 
realistic  and  more  accurate  than  the  Sandia  icsults.  The  Figure  also  verifier  the 
existence  of  a  predominantly  axial  flow  over  most  of  the  cylinder  length,  except 
within  a  region  of  the  order  of  the  radius  near  the  end  wall.  Linear  and  nonlinear 
velocity  distributions  are  hardly  distinguishable.  Clearly  visible  is  the  turning  of 
the  flow  near  the  end  wall.  While  the  flow  appears  steady  in  the  coordinate  sys¬ 
tem  chosen,  the  velocity  field  describes  in  fact  an  oscillatory  motion  of  fluid  ele¬ 
ments  about  their  near-circular  orbit. 

The  pressure  distributions  for  the  same  case  are  shown  in  Figure  3  with  the 
heavy  lines  indicating  positive  values.  Remarkable  is  the  formation  of  regions  of 
high  and  low  pressure  in  the  corner  near  <j>  ~  45°  and  ~  135°  ,  respectively, 
which  produce  large  contributions  to  the  moments  about  i-axis  and  y-axis. 
Except  in  this  region  near  the  end  walls,  the  variation  of  the  pressure  is  relatively 
weak.  The  azimuthal  position  of  extremum  pressure  changes  from  d>  =  0  for 
small  values  of  Re  to  <f>  =  90°  as  Re  — *  oo. 

The  dominant  components  of  velocity  and  pressure  fields  are  azimuthally 
periodic  with  period  2".  The  harmonics  are  small,  indicating  the  small  effect  of 
nonlinearity  in  the  range  of  low  Reynolds  numbers.  The  only  important  non¬ 
linear  term  is  the  aperiodic  mean  flow.  This  is  clearly  shown  by  Figure  1  which 
gives  the  azimuthal  velocity  in  the  center  plane  z  —  0.  The  aperiodic  com¬ 
ponent  is  opposite  to  the  rigid-body  rotation  and  exerts  a  negative  roll  moment 


Appendix  F 


71 


through  the  wall-shear  stress  tt^.  The  axial  and  radial  mean  velocity  field  is 
given  in  Figure  5.  This  streaming  term  exhibits  a  toroidal  motion  near  the  end  in 
each  half  of  the  cylinder  and  causes  a  slow  drift  of  fluid  elements  with  respect  to 
circular  orbits.  This  mean  velocity  produces  the  symmetric  pattern  in  flow  visual¬ 
izations  11  at  low  Reynolds  numbers. 

At  the  higher  Reynolds  number  Re  =  300,  the  maximum  axial  velocity 
appears  at  <t>  ~  90°  .  As  shown  in  Figure  6,  the  flow  in  the  plane  <t>  =  90° 
breaks  up  into  two  swirls,  one  in  each  half  of  the  cylinder,  with  little  flow  across 
the  plane  z  =  0.  Three  weak  swirls  develop  in  the  plane  6  —  0  such  that  the 
velocity  field  is  reminiscent  of  a  chain  with  five  links.  Notably,  the  break-up  into 
cells  is  restricted  to  an  inner  region  of  the  cylinder.  The  motion  in  the  pro¬ 
nounced  boundary  layer  visible  in  the  plane  o  —  0  dees  not  follow  the  cellular 
structure  and  may  have  a  direction  opposite  to  the  core  flow.  The  pressure  varia¬ 
tion  is  characteristically  different  from  that  at  low  Reynolds  number.  Figure  7 
shows  the  strong  variation  and  the  formation  of  an  almost  symmetric  pattern 
along  the  cylinder  in  the  plane  <f>  =  0,  while  the  variation  at  <*)  =  90°  is  rather 
weak.  This  pressure  field  explains  the  void  observations  of  Miller  15  which  show 
a  wavy  distortion  of  the  void  in  the  plane  <p  —  0  at  high  Reynolds  numbers.  The 
free  surface  in  these  observations  can  be  interpreted  as  a  surface  of  constant  total 
pressure.  The  steep  and  opposite  pressure  gradients  across  the  cylinder  axis  near 
z  jr]  —  0.25  and  z  jt]  —  0.75  displace  the  void  near  these  positions  in  opposite 
directions  along  the  diameter  at  <p  ~  -  15°  . 

Calculation  of  the  Liquid  Moments 

Conservation  of  angular  momentum  for  the  steady  flow  in  a  control  volume 
V  with  surface  5  rotating  with  constant  rate  Q  about  a  fixed  axis  requires 

M  +  /  (r  X  F )dS  =  /  r X  (20  Xv)pdV 
s  v 

+  f  r  X  (n  X  (n  X  r)j pdY  +  /  (r  Xv)/;(v VS)  (15) 

V  s 

where  the  velocity  v  is  measured  relative  to  ihe  aeroballistic  frame.  On  the  Icft- 
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hand  side,  M  is  the  resultant  torque  on  the  control  volume,  r  is  the  position  vec¬ 
tor,  and  F  the  stress  acting  on  the  cylinder.  The  presence  and  meaning  of  cer¬ 
tain  terms  depend  on  the  choice  of  the  control  volume.  The  surface  integral  on 
the  right-hand  side  of  Eq.  (15)  vanishes  if  the  surface  of  the  control  volume  is 
closed. 

For  ease  of  practical  application,  we  express  ,  the  moment 
M  =  ( Mx  ,  My  ,  M2 )  in  terms  of  cartesian  components  which  provide  yaw,  pitch, 
and  roll  moment,  respectively.  Analogue  to  Eq.  (1)  we  decompose  the  moments 
into 

M  =  Mr  +  (IB) 

where  Mr  corresponds  to  the  pure  rigid-body  motion  while  originates  from 
the  deviation  velocity  and  pressure.  For  the  cylindrical  control  volume,  the 
rigid-body  rotation  causes  only  a  pitch  component 

m;  -  2*0,11  +  ^-(i  -  |„s)i  (i?) 

while  Ml  —  Ml  —  0.  Note  that  M  is  dimensionless;  the  reference  moment  is 

5  0 

p  a  °oj~.  „ 

The  evaluation  of  the  components  of  bears  some  ambiguity  that  can  be 
exploited  for  advantages.  Previous  computational  work  3>  4>  6  employed  a  con¬ 

trol  volume  consisting  of  an  “empty”  closed  cylinder  with  only  the  pressure  and 
stresses  acting  on  the  inside  surface.  In  this  case,  the  right  hand  side  of  Eq.  (15) 
vanishes  and  the  moments  are  obtained  from  the  stresses  F  at  the  inside  wali  of 
the  control  volume.  Here,  we  use  a  different  choice  that  bears  great  advantages 
especially  for  computational  work. 

We  consider  a  control  volume  consisting  of  a  solid  cylindrical  surface  com¬ 
pletely  enclosing  the  liquid.  The  moment  calculation  for  this  “full”  control 
volume  rests  on  the  relation 

M d  =  J  r  X  (2Q  Xvd)pdV  (18) 

v 

Using  analytical  relations  derived  by  Rosenblat  et  al.,6  the  components  of  Mrf 
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can  be  shown  to  take  the  form 


Mj  =  -  I !  COS# 

CO 

My  =  1 2  sin#  -  1 3  cos# 

(19b) 

Md  =  1 4  sin# 

(19c) 

where 

/j  =  -J4=  j  z(vrcos<f>-  VjSin  <f>)rdrd  <j>  dz 

V 

(20a) 

/2  =  —J  v4)r2drd  <f>dz 

1  V 

(20b) 

J3  =  -  f  vzs\n<f>r2drd  <f>dz 

V 

(20c) 

Finally, 

we  obtain  the  moments  in  the  form 

n  2r  t 

Mx  =  — e—  J  J  f  vzr  2 cos<f>  dr d  <j>  dz 

tan^  -  r)  0  0 

(21a) 

V  2r  1  i)  It  1 

My  —  ej  J  f  v (j>r2drd  <j>dz  H — ~rj  J  f  vzr2s\n<f>drd  <t>dz 

-T)  0  0  tan^  -  rt  0  0 

(22b) 

Mz  =  Mxd  tan# 

(22c) 

The  volume  integral  approach  thus  leads  to  handy  expressions  which  involve  only 
the  radial  and  azimuthal  velocity  components.  Integration  over  d  <f>  reduces  the 
requirements  in  fact  to  the  knowledge  of  the  aperiodic  component  of  v $  and  the 
simply  periodic  components  of  vz .  Therefore,  the  volume  approach  can  also  be 
applied  to  the  analytical  results  given  above  and  provides  yaw  and  pitch 
moments  without  explicit  knowledge  of  the  pressure. 

Results  for  the  Liquid  Moments 

While  velocity  and  pressure  fields  are  primarily  of  basic  fluid  mechanical 
interest,  the  practical  need  for  the  moments  dictates  the  measure  for  efficiency  of 
the  code.  The  moments  derived  from  the  volume  approach  and  the  surface 
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approach  applied  to  the  same  spectral  solutions  are  shown  in  Tables  1  and  2, 
respectively.  The  Reynolds  number  Re  =  20  is  in  the  range  of  maximum  despin 
moment  A f2 . 

It  is  obvious  that  the  volume  approach  provides  results  of  superior  quality 
and  more  rapid  convergence.  The  required  (absolute)  accuracy  of  10" 3  for 
engineering  applications  can  be  achieved  with  the  low  truncation  K  =  4, 
L  =  3,  M  =4.  This  accuracy  has  to  be  seen  in  the  light  of  considerable  uncer¬ 
tainty  in  the  moments  governing  the  exterior  aerodynamics  of  the  projectile.  As 
a  rule  of  thumb,  an  increase  in  the  aspect  ratio  requires  additional  expansion 
functions  in  axial  direction  while  increasing  Reynolds  number  requires  higher 
resolution  in  both  radial  and  axial  direction. 

Figure  8  compares  the  calculated  roll  moments  for  a  wide  range  of  Reynolds 
numbers  with  the  experimental  results  of  Miller  2  and  with  computational 
results.4,  6  The  deviation  of  the  results  of  the  Sandia  code  4  is  caused  by  using 
inappropriate  formulas  for  the  moments  in  the  nutating  coordinate  system.7  The 
agreement  with  the  other  computational  data  is  good.  Test  runs  with  high  reso¬ 
lution  suggest  that  the  small  difference  from  the  results  of  Rosenblat  et  al.6  is  due 
to  lower  resolution  of  the  finite-element  code  in  combination  with  the  application 
of  the  surface  approach  for  the  moments.  The  experiments  were  made  in  a  range 
of  spin  rates  u  between  2000  and  4000  rpm.  While  oj  —  3000  rpm  has  been  used 
in  Figure  8,  assumption  of  a  lower  value  would  improve  the  comparison  with 
respect  to  the  maximum  values. 

Figure  9  shows  a  similar  comparison  for  the  yaw  and  pitch  moments.  The 
results  of  the  Sandia  code  are  suppressed  since  they  suffer  from  a  dimensional 
inconsistency.7  While  the  agreement  for  the  yaw  moment  at  high  Reynolds 
numbers  is  surprisingly  good,  the  deviation  in  the  pitch  moment  is  likely  to  ori¬ 
ginate  from  insufficient  resolution  of  the  steep  pressure  gradients.  This  effect  of 
discretization  errors  has  been  reduced  in  the  spectral  code  by  using  the  volume 
approach  for  calculating  the  moments. 
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Figure  10  show3  the  dependence  of  the  yaw  and  pitch  moments  per  unit 
length  (the  roll  moment  is  proportional  to  Mz  )  on  the  aspect  ratio  of  the  cylinder 
and  compares  with  results  of  the  code  written  by  Strikwerda  Sc  Nagel  7  and  the 
analytical  results  for  r]  — *  oc.  This  diagram  indicates  that  a  reduction  in  the 
overall  liquid  moments  for  a  given  fluid  mass  can  be  achieved  by  splitting  the 
cylindrical  volume  into  slices  or  low  aspect  ratio. 

Discussion 

The  codes  previously  in  use  may  serve  for  establishing  some  basic  results  but 
are  too  inefficient  for  routine  applications.  The  finite-difference  code  developed  at 
Sandia  Laboratories  3|  4  rests  on  Chorin’s  method  of  artificial  compressibility  and 
provides  the  steady  solution  at  11  X  24  X21  grid  points  in  r,  d>,  z -direction  by 
integrating  over  typically  104  to  8  104  time  steps,  a  task  that  requires  6  to  48 
minutes  of  CPU  time  on  a  Cray-lS.  The  result  consists  of  over  22,000  values  of 
the  velocities  vT ,  v  a,  vz  and  the  pressure  p  at  the  grid  points. 

Strikwerda  Sc  Nagel  5  briefly  describe  a  code  using  finite  differences  in  radial 
and  axial  direction  and  pseudospectral  differencing  in  the  azimuthal  direction. 
Nonuniform  grids  are  introduced  for  increased  resolution  near  the  walls.  The 
difference  equations  are  solved  by  an  iterative  method  based  on  successive  over¬ 
relaxation.  The  computer  time  required  is  comparable  to  that  of  the  Sandia 
code.  A  thorough  evaluation  of  the  two  codes  is  currently  conducted  at  BRL.7 

The  experience  with  the  present  version  of  the  spectral  code  shows  that  high 
performance  can  be  achieved.  The  solution  is  obtained  in  serni-analytical  form 
with  only  jY  =  4  K  L  M  (typically  less  than  500)  numerical  coefficients.  This 
low  data  volume  is  especially  attractive  for  storage  and  for  communication  with 
remote  supercomputers.  The  code  is  very  well  suited  for  vcctorizntion.  since 
practically  all  CPU  time  is  spent  on  conslrimting  and  solving  an  algebraic  system. 
The  code  demands  larger  memory  than  other  codes,  because  54-bit  arithmetic  is 
highly  recommended  Cr  spectral  methods  in  general,  and  the  algebraic  system 
requires  A’(,Y  4  1)  words  of  storage.  A  run  with  .Y  =  -  500  requires  aboil*  2 
Mbyte  of  memory  and  can  easily  be  carried  out  on  engineering  workstations 
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within  a  few  minutes,  while  moment  calculations  with  N  =  192  are  a  matter  of 
seconds.  Since  the  memory  requirement  is  acceptable  even  if  higher  resolution  is 
desired,  the  method  applied  here  is  a  viable  alternative  in  numerous  other  fluid 
mechanical  problems.  The  ability  to  obtain  accurate  solutions  for  the  steady 
problem  directly  from  (large)  algebraic  systems  bears  valuable  potential  to  answer 
the  question  whether  the  steady  solution  is  stable,  and  allows  analysis  of 
unsteady  motions  with  implicit  time-stepping.  The  design  of  a  reliable  code  for 
the  unsteady  problem  can  take  profit  from  the  knowledge  of  the  eigenvalue  spec¬ 
trum  for  small  unsteady  disturbances  of  the  stead;  flow. 

While  the  calculation  of  velocity  and  pressure  fields  provides  insight  into  the 
physics  of  the  flow,  the  practical  interest  in  the  moments  for  the  quasi-steadily 
changing  parameters  in  flight  simulations  can  be  satisfied  with  modest  amounts 
of  computer  time.  This  is  due  to  using  a  modified  Newton  method  which  updates 
the  Jacobian  only  when  demanded  by  deteriorating  convergence. 

In  genera!,  the  volume  approach  provides  much  more  accurate  results  than 
the  surface  approach.  This  improvement  is  due  to  the  additional  smoothing  of 
fluctuating  data  by  integrating  over  three  instead  of  two  space  directions  and  to 
us>ng  fewer,  less  fluctuating,  and  more  accurate  input  data.  The  absence  of  vr  in 
the  volume  formulation  is  welcome.  This  velocity  component  is  small  over  most 
ol  the  cylinder  length  but  oscillatory  in  the  radial  direction  10  with  considerable 
gradients  near  the  walk  Near  the  end  walls,  vr  is  of  the  same  order  as  vt  with 
steep  gradients  toward  the  end  wall.  Inspection  of  the  velocity  plots  of  Vaughn 
et  al.4  indicates  that  these  gradients  were  difficult  to  resolve  by  the  finite 
difference  method.  The  aperiodic  component  of  v  $  is  a  relatively  small  streaming 
term  of  smooth  and  almost  uniform  behavior  along  the  cylinder  axis.  The  large 
azimutiially  periodic  components  of  v#  near  the  end  walls  do  not  affect  the 
moment  calculation. 

Probably  the  greatest  advantage  of  the  volume  formulation  is  the  absence  of 
the  pressure  from  the  moment  equations.  This  property  favors  the  use  of 
pressure-free  sets  of  basic  equations,  e.  g.  in  terms  of  vorticity  or  vector  potential. 
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The  smaller  number  of  dependent  variables  can  be  exploited  for  further  increas¬ 
ing  the  efficiency.  Even  in  natural  variable  formulations,  the  pressure  is  difficult 
to  obtain  with  high  accuracy  because  of  the  invalidity  of  the  equations  in  the 
joints  of  the  fiat  end  walls  to  the  cylindrical  side  wall.  As  shown  in  Figure  3,  the 
pressure  may  assume  extrema  near  the  corners  and,  therefore,  inaccuracies  in  this 
legion  may  strongly  influence  yaw  and  pitch  moment.  In  this  context,  it  is 
instructive  to  evaluate  the  convergence  history  of  the  artificially  time-dependent 
method  implemented  in  the  Sandia  code.3  While  the  velocity  rapidly  reaches  a 
quasi-steady  state,  about  75%  of  the  iterations  are  spent  on  improving  the  pres¬ 
sure  field.  We  estimate  that  by  use  of  the  volume  approach  equivalent  or  supe¬ 
rior  values  for  the  moments  could  be  obtained  with  less  than  20%  of  the  itera¬ 
tions.  It  is  worthwhile  to  note  that  the  analytical  results  of  Rosenblat  et  ah, 6 
and  equations  (19),  (20)  for  the  moments  are  valid  for  closed  containers  of  more 
general  shape  and  thus  can  be  used  for  other  interior  flow  problems. 

Our  analytical  and  numerical  tools  allow  quick  estimates  and  efficient  calcu¬ 
lation  of  accurate  liquid  moments.  These  results  also  suggest  guide  lines  for  the 
suppression  of  the  flight  instability  caused  by  the  viscous-liquid  payload.  For  a 
given  cavity  and  fluid,  a  reduction  in  the  overall  liquid  moments  can  be  achieved 
in  two  ways.  The  first  method  is  the  split  of  the  cylindrical  volume  into  slices  of 
lew  aspect  ratio.  The  second  way  is  the  longitudinal  split  into  k~  “straws”  of 
high  aspect  ratio  kq.  The  change  of  the  radius  reduces  the  Reynolds  number 
which  may  or  may  not  be  desirable.  The  nondimensional  moment  per  straw 
increases  due  to  the  increasing  aspect  ratio.  An  essential  reduction  of  the  overall 
moments,  however,  originates  from  the  fact  that  the  dimensional  moments  are 
proportional  to  the  fifth  power  of  the  radius.  The  dimensional  factor  is  therefore 
reduced  by  k~  ‘  per  straw  or  k  3  for  ail  straws  together.  As  a  raw  estimate,  the 
effective  moments  can  be  reduced  by  a  factor  k~ 2. 
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CAPTIONS  OF  TABLES 


Table  1.  Volume  approach  for  the  moments  at  77  =  4.368,  r  =  .1667,  6  =  20, 
Re  ^  20. 

Table  2.  Surface  approach  for  the  moments  at  r?  =  4.368,  r  =  .1667,  6  —  20, 
Re  =  20. 


CAPTIONS  OF  FIGURES 

Figure  1.  Definition  sketch. 

Figure  2.  Vector  plot  of  the  axial  and  radial  velocities  in  the  planes  =  45° 
(left,  scale  0.075)  and  <i>  =  135°  (right,  scale  0.0375)  at  Re  =  20  for  z  >  0. 

Figure  3.  Contour  plot,  of  Die  pressure  field  in  the  planes  <j>  —  45°  (left)  and  and 
4>  —  135°  (right)  at  Re  =  20  for  z  >  0.  Levels  every  0.0025. 

Figure  4.  Vector  plot  of  the  azimuthal  velocity  v $  in  the  center  plane  at  ;  =  0, 
Re  =  20.  Scale  0.003. 

Figure  5.  Vector  plot  of  the  axial  and  radial  mean  velocity  v $  for  Re  —  20. 
Scale  0.002. 

Figure  6.  Vector  plot  of  the  axial  and  radial  velocities  in  the  planes  d  =  0° 
(left,  scale  0.05)  and  0  =  90°  (right,  scale  0.2)  at  Re  =  300  for  c  >  0. 

Figure  7.  Contour  plot  of  the  pressure  field  in  the  planes  0  —  0°  (left)  and 
—  90°  (right)  at  Re  =  300  for  c  >  0.  Le  very  0.005. 

Figure  8.  Roll  moment  M.  >s.  Reynolds  number  Rc  for  rj  =  -1.368, 
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r  =  0.16667,  and  0  —■  20° .  Comparison  with  numerical  and  experimental  data. 


Figure  9.  Yaw  moment  Mz  and  pitch  moment  My  vs.  Reynolds  number  Re  for 
??  —  4.368,  r  =  0.16667,  and  6  —  20° .  Comparison  with  numerical  data. 

Figure  10.  Yaw  moment  Mx  and  pitch  moment  My  per  unit  length  vs.  aspect 
ratio  7?  at  Re  =  10,  r  =  0.16667,  and  6  =  2°.  Comparison  of  present  numeri¬ 
cal  results  with  analytical  results  for  r;  — *  oo  and  data  obtained  by  Nusca  with 

Strikwerda’s  code. 
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Table  1.  Volume  approach  for  the  moments  at  rj  =  4.368,  r  =  .1667,  9  =  20, 
Re  =  20. 
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0.07385 
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Table  2.  Surface  approach  for  the  moments  at  t?  =  4.368,  r  =  .1667,  0  =  20, 
Re  =  20. 
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M, 
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0.07394 
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0.03308 
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0.07247 
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0.07291 

0.03024 
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0.03028 
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0.03036 

6 

5 

8 
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Fig 

(left 


Figure  6.  Vector  plot  of  the  axial  and  radial  velocities  in  the  planes  <$> 
(left,  scale  0.C5)  and  =  90°  (right,  scale  0.2)  at  Re  =  300  for  2  >  0. 
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Figure  8.  Roll  moment  Mt  vs.  Reynolds  number  Re  for  rj  =  4.368 
r  =  0.166C7,  and  0  =  20°  .  Comparison  with  numerical  and  experimental  data. 
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Figure  9,  Yaw  moment  Mz  and  pitch  moment  My  vs.  Reynolds  number  Re  for 
V  =  4.368,  t  =  0.16667,  and  6  =  20°.  Comparison  with  numerical  data. 
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Figure  10.  Yaw  moment  Mx  and  pitch  moment  My  per  unit  length  vs.  aspect 
ratio  t)  at  Re  —  10,  r  —  0.16667,  and  0  =  2°.  Tomparison  of  present  numeri¬ 
cal  results  with  analytical  results  for  r?  — ►  oo  and  data  obtained  by  Nusca  with 
Strikwerda's  code. 
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Abstract 

The  use  of  distributed  workstations  instead  of 
mainframes  has  enabled  us  to  adapt  the  software  to 
our  needs  and  to  implement  symbolic  manipulation 
into  our  research  environment.  Our  applications  are 
primarily  the  simplification  and  verification  of  tedious 
manual  procedures  rather  than  solving  complete  prob¬ 
lems  in  symbolic  form.  We  report  our  cxpencnce 
using  Maesyma  for  the  denvation,  coding,  and  docu¬ 
mentation  of  complex  equations,  development  of 
improved  algorithms  for  spectral  methods,  and  gen¬ 
eration  of  Fortran  code  to  provide  test  data.  Example 
problems  arc  in  the  area  of  viscous  flow  computation 
and  hydrodynamic  stability. 

1.  Introduction 

The  exploration  of  fluid  dynamic  phenomena 
and  their  underlying  mechanisms  often  requires  solv¬ 
ing  the  complex  equations  of  motion  in  relatively  sim¬ 
ple  geometry.  The  arsenal  of  tools  for  solution  ranges 
from  brute -force  computation  to  sophisticated  asymp¬ 
totic  methods  each  of  which  has  merits  and  shortcom¬ 
ings.  Pure  computation  with  finite -difference  or 
finite-element  methods  is  usually  the  shortest  and 
most  direct  route  to  results  but  leaves  the  investigator 
with  vast  amounts  of  data  for  tedious  postprocessing. 
Computation  and  ['cstproccssing  must  he  repeated  for 
every  new  set  of  |  arameters.  Analytical  methods,  cm 
the  other  end.  olicn  provide  generic  results  that  arc 
easy  to  understand  and  clearly  reveal  the  role  of 
parameters.  Derivation  of  analytic.il  results  for  com¬ 
plex  problem',  however,  can  lx:  lime-consuming  and 


prone  to  errors,  and  the  necessary  approximations 
may  blur  the  quantitative  aspects.  Only  in  rare  eases 
can  closed-form  solutions  be  found.  Usually,  analyti¬ 
cal  work  reduces  the  partial  differential  equations  to 
ordinary  differential  equations  which  can  be  solved 
numerically  with  relative  ease. 

In  our  strive  for  quantitative  infomiabon  and 
insight,  we  have  successfully  applied  a  combination 
of  analytical  modeling  and  spectral  methods  for  the 
numerical  work.  Spectral  methods  have  very  attrac¬ 
tive  numerical  properties  1  and  arc  closely  related  to 
the  mathematical  formulation  of  the  problem.  Incon¬ 
sistencies  in  the  approach  do  not  remain  localized, 
e.g.  to  a  tew  grid  points  in  a  finite-difference  solution, 
but  affect  the  overall  solution.  The  sensitivity  of 
spectral  methods  to  even  small  errors  -  in  the  formu¬ 
lation  or  due  to  round-off  -  is  frustrating  for  the 
beginner  but  a  welcome  indicator  for  the  experienced 
user.  Often,  the  inconsistencies  arise  from  inappropri¬ 
ate  treatment  of  singularities,  eg.  comers  of  the 
domain,  from  errors  in  equations  or  code,  or  from 
loss  of  significant  digits  that  may  remain  undetected 
in  other  calculations.  The  successful  effort  of  debug¬ 
ging  the  whole  approach  is  rewarded  with  very 
pleasant  performance  of  the  method. 

The  development  and  implementation  of  correct 
equations  can  he  a  task  of  unexpected  difficulty. 
Repetition  of  derivations  and  coding  is  time- 
consuming  but  not  fail-safe  Hie  a  posteriori  check 
of  results  with  existing  woik  or  intuition  is 
insufficient  for  verification  of  research  compulations. 
We  luxe,  llieiefuie,  incoiporalcd  xymkilic  manipula- 
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tion  with  Macsyma  into  our  studies  on  linear  and  non¬ 
linear  stability  of  shear  flows  and  on  viscous  rotating 
flows. 

In  our  applications,  Macsyma  has  been  used  less 
for  problem  solving  than  for  reproducing  and  verify¬ 
ing  the  steps  of  the  analysis,  Like  most  computer 
tools,  Macsyma  requires  major  learning  efforts,  while 
tutorial  examples  for  a  specific  step  in  a  practical 
application  arc  hard  to  find.  With  increasing  skills, 
however,  the  efficiency  in  performing  symbolic  mani¬ 
pulation  dramatically  increases.  The  gain  in  overall 
efficiency  is  high  enough  to  consider  new  classes  of 
problems  and  new  approaches  that  would  hardly  be 
feasible  without  symbolic  manipulation.  We  have 
also  experienced  that  this  tool  trains  the  user  to 
clearly  formulate  the  steps  of  the  analysis  and  to 
"play”  with  alternatives  which  at  times  provide  new 
insight  into  the  structure  of  the  problem  and  stimulate 
more  elegant  formulations. 

1.  The  Computational  Environment 

We  have  Macsyma  Version  309.6  installed  on  a 
Sun  3/1 80S  with  575MB  disk.  This  server  is  con¬ 
nected  to  Ethernet  and  to  the  campus-wide  Sonnet 
fiber-optic  link.  Therefore,  we  can  access  Macsyma 
locally  and  remotely  from  other  workstations  or  dial 
in  through  the  university’s  Micom  lines  from  home. 
The  system  is  currently  used  by  about  ten  faculty  and 
graduate  students  with  little  or  no  prior  experience  in 
Lisp  or  any  symbolic  manipulators.  Some  members 
of  this  group  have  worked  with  Unix-based  worksta¬ 
tions  and  superminis  for  the  past  five  years  and  had  a 
short  encounter  with  Macsyma  on  a  Masscomp 
workstation  with  insufficient  disk  space.  The 
Macsyma  software  is  not  under  maintenance,  there¬ 
fore,  tire  reference  manual  2  and  users  guide  3  are  the 
only  common  sources  of  information  for  self-study. 
Although  our  applications  vary  over  a  wide  range  of 
topics  in  fluid  mechanics  and  heat  transfer,  we  have 
not  yet  fully  explored  the  capabilities  of  Macsyma. 

Our  primary  applications  of  symbolic  manipula¬ 
tion  are  in  the  following  areas: 

( 1 )  Derivation  of  equations 

(2)  Generation  of  Fortran  code 

(3)  Comparison  and  selection  of  algorithms 

(4)  Documentation 

We  outline  iri  the  following  sections  the  major  steps 
involved  in  solving  some  sample  problems  and  indi¬ 


cate  areas  of  application  for  symbolic  manipulation, 

2.  The  Flow  in  a  Spinning  and  Nutating  Cylinder 

Calculation  of  the  flow  in-  a  spinning  and  nutat¬ 
ing  cylinder  is  a  prototype  problem  arising  in  the 
design  of  spin-stabilized  rockets  or  shells  with  liquid 
payloads.  Interest  is  in  the  moments  exerted  by  the 
internal  fluid  motion  and  their  detrimental  effect  on 
flight  stability.  Since  the  angle  between  spin  axis  and 
trajectory  is  small  and  the  rate  of  nutation  is  smaller 
than  the  spin  rate,  a  small  dimensionless  parameter  e 
can  be  identified  that  governs  the  deviation  of  the 
fluid  motion  from  rigid-body  rotation.4  This  viscous 
secondary  flow  in  the  cylinder  and  the  coupling  with 
inviscid  inertial  waves  arc  analyzed  using  a  perturba¬ 
tion  method  for  long  cylinders,  an  expansion  in  spa¬ 
tial  eigenfunctions,  and  a  Navicr-Stokes  solver  based 
on  a  spectral  collocation  method.  All  these  methods 
rest  on  the  continuity  equations  and  momentum  equa¬ 
tions  for  the  reduced  pressure  and  the  secondary  velo¬ 
cities  (deviation  from  rigid-body  motion).  The  equa¬ 
tions  are  written  in  cylindrical  coordinates  with 
respect  to  a  rotating  coordinate  system  4  In  this  sys¬ 
tem,  the  secondary  motion  appears  steady  arid  forced 
by  an  azimuthally  periodic  force  of  order  0(E). 

This  set  of  equations,  the  characteristics  of  the 
variables  ( scalar ,  real),  dependencies,  and  frequently 
occuring  operators  can  be  defined  once  and  saved  in  a 
file.  The  branching  into  different  approaches  suggests 
a  hierarchical  file  structure  as  in  Unix.  We  found  it 
advantageous  -  especially  in  the  learning  phase  -  to 
develop  segments  of  this  hierarchy  interactively  and 
to  save  the  successful  commands  in  relatively  small 
files  which  terminate  with  cleanup  commands  (kill). 
This  mode  is  very  conveniently  executed  on  multi- 
window  systems  with  easy  scroll  and  scrccn-to-screcjt 
copy.  The  proper  sequence  of  files  to  be  loaded  can 
be  specified  in  a  master  file. 

2.1  Perturbation  Analysis  for  Long  Cylinders 

For  a  sufficiently  long  cylinder,  the  soluiipn  can 
be  assumed  independent  of  the  axial  direction 'and  the 
boundary  conditions  on  the  end  walls  can  be 
neglected.  The  essential  steps  of  the  perturbation 
expansion  arc:  specify  axial  gradients  to  be  zero  (gra¬ 
de/),  expand  velocities  and  pressure  in  a  truncated 
power  series  in  e  (sum),  expand  azimuthally  in  a  com¬ 
plex  Fourier  series,  and  substitute  the  double  sum  into 
the  basic  equations.  Extraction  of  tcrnis  in  the  nth 
power  of  c  and  like  exponentials  e'"’0  in  the  aximu- 
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thal  variable  $  with  raicoeff(e\(...))  provides  ordinary 
differential  equations  for  the  coefficient  functions 
which  depend  on  the  radial  variable  r.  This  pro 
ccdurc  shows  an  interesting  pattern,  since  numerous 
equations  arc  identically  satisfied:  odd  n  arc  associ¬ 
ated  with  odd  m  and  contribute  only  to  the  axial  velo¬ 
city  component,  while  even  n  arc  associated  with 
even  m  and  contribute  to  radial  and  azimuthal  veloci¬ 
ties  and  the  pressure.  Wiii;  homogeneous  boundary 
conditions  on  the  cylinder  wall  and  proper  conditions 
on  the  axis,  the  lowest-order  axial  velocity  component 
can  be  obtained  in  closed  form  and  expressed  in  terms 
of  the  modified  Bessel  function  l \(<xr),  where 
a"  -  iRc  and  Re  is  the  Reynolds  number.  The  mean 
velocity  {m  even)  is  of  order  O  (c:)  and  consists  only 
of  an  azimuthal  component.  This  component  is  asso¬ 
ciated  with  a  viscous  roll  moment  on  the  cylinder. 

We  have  manually  developed  a  closcd-lorm 
solution  in  terms  of  Bessel  functions  for  the  periodic 
second-order  terms.5  We  have  not  yet  reproduced 
this  result  with  Macsyma.  The  effort  to  obtain 
closed-form  solutions  increases  dramatically  with  the 
order.  Therefore,  we  solve  for  the  terms  of  order 
n  >  2  by  use  of  a  spectral  collocation  method.  The 
expansion  functions  arc  linear  combinations  of  even 
and  odd  Chcbyshcv  polynomials  in  the  interval 
0  <  r  <  1  that  satisfy  the  boundary  conditions  3nd 
conditions  at  r  =0.  The  case  n  =  2  serves  for  com¬ 
parison  with  the  closed- form  solutions. 

Given  the  (linear)  differential  equation  Lf  -  g , 
where  /  -f(x)  and  g  -  g(x)  and  the  expansion 
functions  hk(x),  spectral  collocation  directly  leads  to  a 
linear  algebraic  system  Ma  =  g.1  The  vector  g  con¬ 
sists  of  the  values  g  (xy )  at  the  collocation  points  x;. 
In  our  work,  the  vector  a  consists  of  the  expansion 
coefficients  ak  of  /  with  respect  to  hk(x),  and 
M  =  (mjk )  is  the  matrix  representation  of  the 
differential  operator  applied  to  hk  (columns)  at  xj 
trows).  The  key  to  solving  the  given  differential 
equation  is  to  accurately  define  the  elements  mjk  of 
die  matrix  M  and  the  right-hand  side  g  The  alge¬ 
braic  system  can  then  be  solved  with  standard  pro¬ 
cedures. 

With  Macsyma,  it  is  straightforward  to  replace 
the  right  hand  side  g  by  a  vector  form  g \n )  and  to  do 
llic  same  for  the  r  dependent  coefficients  of  the 
differential  equation  What  remains  then,  is  to  replace 
iltff  if  .r  ,n)  by  ll'j.k.n).  n  >0,  where  j  indicates 
the  collocation  fxiini,  k  die  expansion  function,  and  n 


the  order  of  differentiation.  The  result  can  be  con¬ 
verted  to  Fortran  and  inserted  into  a  generic  Fortran 
framework.  The  Macsyma  output  usually  requires 
some  editing  for  inconsistencies  (e.g.  sqrt(2 ))  with  our 
compiler.  After  this  editing  and  modifications  to 
account  for  the  specific  parameters  of  the  problem, 
printout,  etc.,  the  spectral  code  is  ready  to  run,  pro¬ 
vided  the  array  H(j,k,n)  contains  the  required  data. 
These  data  are  generated  by  one  of  a  number  of 
problem-independent  subroutines  in  our  individual 
library.  Different  versions  account  for  different  sym¬ 
metries,  intervals,  variable  transformations,  and  the 
specific  choice  of  collocation  points  (see  section  4). 

2.2  Volume  Integration  of  the  Moments 

Given  the  deviation  velocities  and  reduced  pres¬ 
sure,  the  calculation  of  the  moments  exerted  by  the 
fluid  on  the  cylinder  is  a  formidable  task.  In  previous 
finite-difference  codes,  these  moments  arc  obtained  by 
integrating  normal  and  tangential  stress  components 
over  the  surface  of  the  cylinder.  Unfortunately,  the 
input  data  -  surface  pressure  and  velocity  gradients  - 
arc  not  very  accurate.  In  addition,  this  procedure 
obtains  the  moments  as  the  small  difference  of  large 
numbers  (pressure  and  viscous  contribution)  and  con¬ 
sequently  the  results  suffer  from  insufficient  resolu¬ 
tion. 

With  some  analytical  effort,  the  moments  can  be 
obtained  by  integrating  over  the  volume  of  the 
cylinder.  The  final  equations  given  by  Herbert  &  Li  6 
show  that  neither  the  pressure  nor  derivatives  of  the 
velocity  arc  required;  only  the  mean  azimuthal  velo¬ 
city  and  the  first  Fourier  mode  (-e"5)  of  the  axial 
vcloc  iy  need  to  be  known.  The  expressions  do  not 
involve  small  differences  of  comparable  terms. 

The  original  derivations  of  the  moment  equa¬ 
tions  by  Rihua  Li  consumed  weeks  of  intense  work. 
The  smaller  part  of  this  time,  obviously,  served  to 
find  the  proper  steps  of  the  analysis.  The  remainder 
was  used  to  derive  and  verify  the  detailed  equations 
and  to  insert  and  correctly  integrate  the  spectral 
representation  of  the  velocity  components.  Since 
Chcbyshcv  polynomials  arc  closely  related  to  powers 
of  the  variable  x  as  well  as  to  trigonometric  functions 
of  the  transformed  variable  0  =  cos_lx,  the  spectral 
moment  equations  give  the  result  in  an  casy-lo- 
intcrprcl  semi-analytical  form.  Meanwhile,  we  have 
repealed  I.i's  derivation  with  Macsyma  in  a  lew 
hours. 
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2.3  Expansion  in  Spatial  Eigenfunctions 

For  cylinders  with  smaller  aspect  ratios,  above 
pcr.urbation  solution  deteriorates.  The  problem  is 
then  to  solve  the  basic  equations  for  a  finite  cylinder. 
While  previous  finite-difference  codes  suffer  from 
insufficient  resolution  and  lack  of  convergence  at 
higher  Reynolds  numbers,  the  problem  appears  to  be 
well  suited  for  a  spectra!  approach.  Two  alternatives 
to  consider  arc:  (1)  spectral  expansions  in  generic 
functions,  with  azimuthal  Fourier  series,  radial  and 
axial  series  in  combinations  of  Chcbvshcv  polynomi¬ 
als  (even  or  odd  to  satisfy  the  conditions  at  r  =  0  and 
to  exploit  the  symmetries  of  the  problem),  or  (2) 
expansion  in  azimuthal  Founcr  senes  and  spatial 
eigenfunctions  that  can  be  obtained  from  the  linear¬ 
ized  problem. 

Expansions  in  spatial  eigenfunctions  to  solve  the 
linear  problem  were  pursued  by  Hall  ct  al.7  who 
determined  the  eigenfunctions  in  a  separate  step  by 
numerical  intcgiation  of  homogeneous  ordinary 
differential  equations  in  the  axial  variable  z.  The 
three-dimensional  boundary  value  problem  reduces  in 
this  way  to  the  one-dimensional  problem  of  satisfying 
the  boundary  conditions  at  the  end  wall  with  a  collo¬ 
cation  or  least  square  method. 

l.i  &  Herbert  8  have  developed  a  single  partial 
differential  equation  for  the  spatial  eigenfunctions  that 
can  be  solved  by  separation  of  variables.  The  general 
solution  consists  of  the  product  of  cosines  of  complex 
arguments  in  the  axial  and  modified  Bessel  funclioas 
of  complex  arguments  in  the  radial  direction.  Eigen- 
solutions  can  be  found  by  determining  the  constant  of 
separation  from  a  characteristic  system  of  equations. 
The  remaining  one-dimensional  problem  here  is  to 
satisfy  tlnec  rattier  complicated  boundary  conditions 
at  the  cylinder  wall.  With  the  numerical  expansion 
coefficients  determined  by  collocation  or  least  square 
method,  the  moments  of  tiic  linear  approach  are 
obtained  in  semi-analytical  form. 

The  Holes  and  derivations  lor  this  approach  fill  a 
file  nl  well  nvei  200  pages.  The  results  arc  still 
uri| -til  >ii  si  nil  sime  a  sm.V  error  was  detected  in  the 
final  results  for  the  moments.  Redcrivation  of  the 
complete  formulation  with  Maesyma  by  a  relatively 
unexperienced  sturiciii  required  two  days  and  showed 
one  sign  com  in  the  manual  derivation  that  led  to 
cancel!. iliori,  not  summation  of  two  small  terms  At 
the  end  ol  the  second  day,  we  obtained  die  correct 
result  hum  the  |  nrn.wi  cxlc. 
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The  expansion  in  spatial  eigenfunctions  is  very 
interesting  from  a  fluid  dynamic  point  of  view  since  it 
clearly  reveals  inertial-wave  type  functions  and  the 
formation  of  boundary  layers  at  higher  Reynolds 
numbers.  However,  the  evaluations  of  convergence 
and  computational  demand  lead  to  the  conclusion  that 
spatial  eigenfunction  expansions  may  not  be  the 
optimum  approach.  Reliable  calculation  of  the  roots 
of  the  characteristic  equations  in  the  large  range  of 
Reynolds  numbers  (10'1  <,  Re  <  106,  say)  is  numeri¬ 
cally  very  demanding  and  can  only  be  achieved  after 
detailed  studies  on  the  asymptotic  form  of  the  eigen¬ 
value  spcctmm.  Symbolic  manipulation  of  the  equa 
lions  and  the  ultimate  solve  command  simplified  this 
task  enormously.  Currently,  we  obtain  reliable  and 
accurate  results  up  to  Reynolds  numbers  of  the  order 
103.  While  we  have  mastered  the  calculation  of  the 
eigenvalues  up  to  Re  -  5  1 06,  more  analytical  refor¬ 
mulation  is  required  to  avoid  floating  point  overflows 
and  loss  of  significant  digits  in  the  numerical  solution 
for  the  expansion  coefficients.  Based  on  existing 
batch  files,  this  work  will  be  done  in  an  interactive 
mode  with  Maesyma  alternating  with  sketching  new 
attempts  on  the  note  pad. 

2.4  Spectral  Navier-Stokes  Solution 

From  previous  work  we  have  learned  that  the 
solution  of  the  cylinder  problem  is  essentially 
governed  by  linearized  equations  plus  nonlinear 
corrections  that  increase  with  e.  Linearization  can  be 
performed  in  different  ways.  The  first  is  a  lineariza¬ 
tion  in  e.  A  second  linear  system  can  be  obtained  by 
linearization  of  the  velocity  components  and  pressure 
about  some  known  solution  which  may  be  identically 
zero.  The  latter  procedure  is  very  efficient  if  the  solu¬ 
tion  is  sought  for  a  densely  spaced  sequence  of 
parameter  combinations  as  in  (light  simulations. 

The  algebraic  form  of  the  equations  is  obtained 
by  use  of  spectral  collocation  with  generic  functions.6 
Instead  of  a  single  array  H(j,k,n)  above,  various 
preset  arrays  //„  //,$,  //„  arc  required,  where  i  indi¬ 
cates  the  dependent  variables.  The  spcctm!  form  of 
the  basic  equations  is  not  difficult  to  codi  manually 
but  -  with  the  Maesyma  input  available  -  can  be 
quickly  generated  per  computer.  The  generic  frame¬ 
work  reflects  the  increased  number  of  dimensions  and 
variables  as  well  as  the  nonlinear  nature  of  the  pint  - 
1cm:  the  dimension  of  the  matrix  increases  to  the  pro¬ 
duct  of  the  numbers  of  expansion  functions  in  r,  <f>,  z 
times  lout  (variables)  The  linear  algebraic  system 
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takes  ihc  form  Mva|A  =  that  is  obtained  by  applying 
Newton’s  method  (p  =  v+1)  or  the  modified 
Newton's  method  (p.>v+l).  In  the  latter  case,  the 
Jacobian  Nlv  is  only  updated  when  indicated  by 
deteriorating  rate  of  convergence.  The  solution  pro¬ 
cedure  starts  wiih  v  =  0,  |i  =  1  from  zero  or  from  a 
previously  obtained  solution  at  neighboring  parame¬ 
ters.  Once  the  iterative  solution  converges,  the 
moments  arc  calculated  from  volume  integrals. 

In  the  Reynolds  number  range  covered  by  previ¬ 
ous  finitc-diflcnence  codes,  accurate  moments  can  be 
obtained  with  4x3x4  expansion  functions  in 
r,  z,  respectively,  and  a  matrix  dimension  of  192. 
The  measured  computer  times  arc  smaller  than  10“3 
times  those  required  by  finite-difference  codes.  These 
threc-ordcrs-of-magniiudc  savings  in  computing  lime 
and  cost  permit  routine  application  of  the  spectral 
code  in  (light  simulations  performed  for  design  pur¬ 
pose  which  was  not  possible  with  previous  codes.  The 
cost  per  solution  increases  with  the  Reynolds  number 
but  is  still  less  than  with  eigenfunction  expansions  at 
Re  ~  103  and  is  for  the  nonlinear  result,  not  the  linear 
approximation.  For  higher  Reynolds  numbers, 
methods  based  on  the  boundary-layer  approximation 
are  effective  although  their  accuracy  at  Reynolds 
numbers  of  the  order  103  remains  to  be  verified. 

The  overall  design  of  the  spectral  code  and  the 
volume  approach  to  the  calculation  of  moments  was 
made  possible  only  by  good  understanding  and 
intense  exploitation  of  the  analytical  properties  of  the 
problem.  Unavailability  of  any  symbolic  manipulator 
in  our  earlier  computer  environment  had  forced  us  to 
do  most  of  the  work  manually.  In  retrospect,  sym¬ 
bolic  manipulation  had  easily  saved  one  year  of 
qualified  labor  on  this  problem  alone.  We  have  pur¬ 
chased  the  license  for  Maesyma  for  $1,800  (for  edu¬ 
cational  institutions),  to  be  shared  by  various  users  -  a 
worthwhile  investment. 

3.  Hydrodynamic  Stability 

The  situation  in  our  studies  of  hydrodynamic 
instability  and  transition  in  shear  flows  is  quite  similar 
to  the  above.  Decomposition  of  the  velocity  and 
pressure  fields  into  various  com|xincnts,  Fourier 
decomposition  in  two  spatial  directions,  and  extraction 
of  linearized  or  nonlinear  equations  for  single  com- 
[xincnis  are  tedious  tasks.  Even  more  tedious  arc  the 
common  steps  of  eliminating  tiic  pressure  by  taking 
the  curl  ol  die  momentum  equations  and  eliminating 
one  ol  ihc  velocity  components  by  use  of  continuity. 


These  steps  require  manipulations  of  derivatives  of 
the  basic  equations  and  subsequent  substitutions.  We 
.me  made  much  progress  in  performing  these  steps 
with  Maesyma  although  elegant  commands  for  some 
detailed  steps  have  not  yet  been  found.  In  our  pro¬ 
cedures,  we  prevent  Maesyma  from  evaluating 
differentiations  to  maintain  differential  equations  in 
the  original  variables  (not  their  derivatives).  Since 
Maesyma  distinguishes,  e.g.  between  (ux)y  and  ( uy)x , 
a  formal  substitution  for  ux  and  subsequent  evaluation 
ev(%,dijf)  does  not  always  provide  the  result  we 
desire.  Our  skills,  though,  arc  still  open  to  improve¬ 
ment. 

A  particular  example  for  the  fruitful  application 
of  symbolic  manipulation  is  our  study  of  nonlinear 
secondary  instability  in  boundary  layers.9,10  In  spite 
of  thorough  checks  of  equations  and  Fortran  code,  we 
ended  up  with  two  independent  codes  that  provided 
different  results.  In  view  of  some  ambiguity  in  cer¬ 
tain  substitutions  and  arrangement  of  terms,  we  were 
unable  to  find  the  source  of  the  discrepancy  and  to 
decide  for  either  one  of  the  codes  to  be  correct. 
While  we  were  caught  for  months  in  this  dilemma 
and  faced  a  third  and  fourth  pass  through  the  analysis, 
Maesyma  became  available  for  MC680OO-bascd 
workstations.  Although  we  had  to  install  and  run  the 
software  on  a  workstation  with  insufficient  disk  space 
and  had  to  learn  the  use  of  Maesyma  from  the  very 
basics,  the  correct  code  was  obtained  within  a  few 
weeks.  The  printout  of  the  equations  and  detailed 
output  from  this  code  enabled  correction  of  the  manu¬ 
ally  derived  equations  and  of  the  faster  hand-written 
code  shortly  after. 

Weaponed  with  the  new  capabilities,  we 
currently  implement  non-parallel  effects  and  the 
streamwise  variations  of  the  disturbances  into  the 
analysis.  Besides  the  creation  of  the  matrices  and 
vectors  for  the  algebraic  systems,  symbolic  manipula¬ 
tion  is  used  to  analyze  asymptotic  properties  of  the 
solulioas,  e  g.  at  large  distance  from  the  wall.  Out¬ 
side  the  boundary  layer,  the  coeflicicnts  of  the 
differential  equations  arc  coastant.  Ihc  asymptotic 
solution  can,  therefore,  be  obtained  in  closed  form. 
Recently,  we  discovered  that  the  existing  theories  of 
stability  in  non-parallel  Hows  use  an  incorrect  form  of 
these  asymptotic  solutions. 

The  complexity  of  the  equations  and  difficulties 
of  the  code  verification  markedly  increase  when  we 
change  from  the  previous  incompressible  problems  to 
similar  problems  in  supersonic  and  hypersonic  flows. 
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The  main  reasons  are  the  increased  number  of  vari¬ 
ables  and  functional  relations,  e.g.  between  viscosity 
and  temperature  (or  species  distribution).  In  joint 
work  with  Dr.  G.  Erlebachcr  at  NASA  Langley,  we 
[vrfonned  the  first  critical  step:  to  develop  and  verify 
a  spectral  code  for  compressible  stability  studies.  The 
verification  used  Maesyma  to  generate  the  matrix  ele¬ 
ments  ol  the  high-dimensional  linear  eigenvalue  prob¬ 
lem.  With  little  effort,  the  effect  of  approximations 
ciiip!o>cd  in  other  codes  could  be  studied  and  a  small 
eiTor  in  another  available  code  discovered. 

Over  time,  the  use  of  symbolic  manipulation 
pi  educes  a  variety  of  relatively  small  input  files  that 
arc  wonhwhile  saving  and  developing  into  an  archive 
of  procedures  comparable  to  a  subroutine  library. 
Improved  versions  can  be  created  as  skills  develop. 
New  projects  can  lake  profit  of  this  archive  of  rcoc- 
cunng  sets  of  equations  or  formal  manipulations  mak¬ 
ing  the  application  of  Maesyma  increasingly  efficient. 

4.  Algorithms  for  Spectral  Methods 

While  Lite  foregoing  applications  were  primarily 
concerned  with  the  derivation  of  error-free  equations 
and  codes,  there  is  a  wealth  of  algorithms  dial  arc 
mathematically  correct  but  useless  for  numerical 
apnlications  -  often  because  of  the  finite  word  length 
foi  handling  numbers  in  a  computer.  The  finite  word 
length  and  round-off  errors  arc  of  particular  concer. 
in  context  with  spectral  methods  and  their  sensitivity 
to  small  errors.  It  is  inevitable  to  secure  accurate 
numerical  results,  especially  for  those  subroutines  that 
produce  the  problem-independent  spectral  data  such  as 
the  matrices  H(j,k,n )  mentioned  above. 

An  early  example  of  a  numerical  pitfall  with 
spectral  methods  is  associated  with  the  first  spectral 
solution  of  the  Orr-Sommcrfcld  equation  of  the  hydro- 
dynamic  stability  theory  by  Orszag.11  This  woik 
employed  the  tau  method  instead  of  collocation  but 
rests  on  matrices  similar  to  H(J,k.n)  to  express  the 
derivatives  of  the  expansion  functions  in  terms  of 
expansion  functions.  For  the  fourth  derivative,  the 
matt ix  elements  lake  the  form 

k[k\k2  -  4)2  -  3 j2kA  +  3jAk2-j2(j 2  -  4)2] 

l  or  j  -  k ,  the  result  of  evaluating  this  expression 
with  32-bit  floating  point  arithmetic  deteriorates 
rapidly  for  j ,  k  >  25,  as  shown  by  Orszag's  results 
for  the  eigenvalues  of  the  Orr-Sommcrfcld  equation. 
The  obvious  reason  is  die  small  difference  between 
the  large  iiutnlx  is  k('  and  2’.  The  effect  of  round-off 


can  be  easily  studied  with  Maesyma  by  obtaining  first 
the  accurate  result  and  then  the  result  of  biguoat 
operations  with  different  sellings  of  fpprec,  the 
desired  precision.  It  is  also  easy  to  see  that  a  numeri¬ 
cally  more  advantageous  form  of  the  matrix  elements 
can  be  obtained  by  substituting  k  -  I  +  j: 

(/  +  j)[l\l2-%)  +  jl\6l2-32) 

+  (12 j2l2  +  8 j3/)(/4  -  4)  +  16/2  +  32y'/) 

This  substitution  is  a  trivial  symbolic  step  but  the 
effect  of  round-off  errors  is  practically  removed. 

Every  step  in  the  subroutines  for  initializing  the 
matrices  H(J,k,n)  has  been  carefully  selected  to 
minimize  round-off  errors.  Calculations  of  data  for  the 
sequence  of  trigonometric  functions  or  Chebyshev 
polynomials  are  often  recursive  and  bear  the  danger 
of  dramatic  error  propagation.  Error  analyses  in  the 
literature  and  comparison  of  results  for  different 
floating-point  precision  were  earlier  used  as  criteria. 
More  recently,  these  studies  have  been  repeated  with 
Maesyma,  although  no  improvements  were  achieved. 
Accuracy  of  the  basic  subroutines  has  been  verified 
for  up  to  180  expansion  functions.  For  larger 
numbers  of  functions  (which  are  irrelevant  for  our 
purpose),  round-off  errors  slowly  creep  in. 

5.  Documentation 

Availability  of  workstations  and  virtually  unlim¬ 
ited  access  to  a  computer  (and  laser  printer)  at  no 
nominal  charge  has  modified  our  mode  of  operation 
over  the  past  few  years.  Besides  program  develop¬ 
ment  and  compulation,  the  computer  also  serves  for 
documentation  and  as  a  “note  book.”  While  docu¬ 
mentation  earlier  involved  piles  of  scratch  notes,  pro¬ 
gram  listings,  and  results  before  the  first  draft  of  a 
publication,  we  usually  begin  this  draft  before  the  first 
line  of  coding.  The  use  of  Unix  has  led  us  to  still 
adhere  to  trofj  as  the  common  typesetting  system 
(although  TeX  is  available).  With  the  typese::  true 
command,  intermediate  steps  and  results  of  our 
Maesyma  scripts  can  be  directly  converted  to  eqn  I 
troff  form  for  inclusion  in  the  descriptive  text.  For 
lengthy  expressions,  die  Maesyma  output  frequently 
requires  some  manual  modification.  The  short  nota¬ 
tion  of  the  scripts  also  suggests  global  redefinition  of 
the  quantities  for  trvff  output.  The  printout  of  the 
Maesyma  batch  script  provides  all  the  details  of  the 
derivations  The  draft  can  be  extended  as  the  work 
progresses  and  provides  an  unambiguous  description 
of  the  equations  used  in  the  computer  codes  Conver- 
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sion  to  a  publication  requires  major  editing,  more  text, 
fewer  equations,  inclusion  of  the  results,  etc.  The 
procedure  prevents,  however,  the  discrepancies 
between  what  has  been  done  and  what  has  been  pub¬ 
lished.  Symbolic  manipulation  serves  here  to  save 
tedious  typing  of  equations  and  to  reduce  sources  of 
error. 

6.  Conclusions 

Symbolic  manipulation  does  not  solve  our 
research  problems.  We  still  have  to  develop  the  ideas 
and  work  out  the  details  of  the  analysis.  However, 
the  description  of  the  various  steps  can  be  sketchy, 
arid  the  accurate  formulation  for  documentation  and 
computer  code  can  be  obtained  by  symbolic  manipu¬ 
lation.  Mnesyma  can  also  be  used  as  a  powerful  and 
efficient  "calculator”  to  test  and  improve  numerical 
algorithms. 

We  consider  symbolic  manipulation  as  a  tool 
similar  to  programming  languages  and  computers:  ini¬ 
tially,  their  use  helped  to  reduce  the  amount  of  human 
labor  and  human  errors.  Soon,  these  tools  enabled 
solutions  to  problems  that  were  not  feasible  without 
their  use.  Learning  and  mastering  these  tools  require 
some  clfon,  but  their  use  educates  us  to  more 
rigorously  formulate  ideas  and  relieves  us  from  time- 
consuming  tasks.  Today,  symbolic  manipulation 
should  be  an  integral  pait  of  the  engineer’s  research 
environment. 
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ABSTRACT 

The  moments  exerted  by  the  fluid  motion  in  liquid  payloads  can  jeopardize  the  stable  flight  of  spin- 
stabilized  projectiles.  These  moments  can  be  computed  with  various  Navier-Stokes  codes  for 
sufficiently  small  Reynolds  numbers.  For  the  linearized  problem.  Hall,  Sedney,  and  Gerber  (1987) 
have  proposed  an  expansion  in  radial  eigenfunctions  that  could  extend  the  range  of  Reynolds 
numbers.  The  need  to  determine  eigenvalues  and  eigenfunctions  numerically  and  the  slow  conver¬ 
gence  of  the  expansion  series,  however,  make  this  approach  less  efficient  than  solving  the  Navier- 
Stokes  equations  with  spectral  methods.  We  have  derived  a  new  formulation  of  the  linear  problem 
that  permits  closed-form  solutions  for  radial  or  axial  eigenfunctions  and  calculation  of  the  eigenvalues 
from  a  closed-form  characteristic  equation.  The  formulation  is  also  suited  for  solution  by  spectral 
methods  at  a  fraction  of  the  computational  expense  of  other  codes. 


L  Introduction 

Spin-stabilized  projectiles  with  liquid  payloads  may  experience  severe  flight  instabilities.  Two 
types  of  liquid-induced  instability  arc  currently  known.  Both  types  are  excited  by  the  coning  motion 
of  the  projectile  about  the  trajectory.1  The  first  type  is  caused  by  resonance  with  menial  waves  at 
critical  frequencies  x  (ratios  of  coning  rate  Q  to  spin  rate  oo).  This  resonance  is  most  pronounced  for 
low-viscosity  liquid  fills,  i.e.  at  high  Reynolds  numbers,  and  depends  sensitively  on  the  cylinder’s 
aspect  ratio  p.  Theoretical  analysis  of  this  instability  usually  involves  the  boundary-layer  approxima¬ 
tion  and  therefore  provides  design  criteria  for  sufficiently  large  Reynolds  numbers,  Re  >  1000.  say. 
Analysis  based  on  the  Navier-Stokes  equations  2  shows,  however,  that  resonance  with  inertial  waves 
may  severely  influence  the  liquid  moments  at  Reynolds  numbers  as  low  as  Re  =  100.  (We  define  the 
Reynolds  number  by  Re  =  pcoa‘/p,  where  p  is  the  density,  a  the  radius  of  the  cylindrical  cavity,  and 
4  the  viscosity.)  The  second  type  of  instability  is  of  viscous  nature  and  is  most  pronounced  at  low 
ar.d  medium  Reynolds  numbers  for  a  wide  range  of  aspect  ratios  and  frequencies.  This  instability  is 
characterized  by  a  negative  roll  moment  that  opposes  the  spinning  motion  of  the  projectile. 
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Theoretical  analysis  of  the  moments  must  be  based  on  the  Navier-Stokes  equations  or  on  a  linearized 
form  of  these  equations  if  the  nutation  angles  9  are  sufficiently  small. 

Since  both  types  of  instability  appear  simultaneously  for  a  wide  range  of  Reynolds  numbers,  it 
would  be  highly  desirable  to  develop  a  unified  approach  that  bears  potential  for  computing  the 
moments  in  a  wide  range.  Moreover,  the  efficiency  of  the  moment  computation  is  crucial  since  flight 
simulations  requite  either  large  tables  for  interpolation  or  frequent  evaluation  of  the  moments  for 
varying  parameters.3 

The  analysis  of  liquid  moments  for  engineering  design  is  based  throughout  on  the  quasi-steady 
motion  in  the  aeroballistic  system.  Direct  solution  of  the  three-dimensional  nonlinear  Navier-Stokes 
problem  typically  requires  5  -  200  seconds  per  solution  on  today’s  workstations  (Sun  3/140  FPA) 
when  using  spectral  methods.2  The  time  increases  with  the  Reynolds  number  and  the  formation  of 
boundary  layers.  For  practical  purpose,  the  application  ranges  up  to  Re  =  1000.  For  larger  Reynolds 
numbers,  Re  >  1000,  the  analysis  rests  on  linearized  equations  and  use  of  the  boundary-layer  approx¬ 
imation4  At  Re  =  OflO3),  this  approximation  causes  errors  of  about  10%  5  which  decrease  with 
increasing  Re.  Estimates  of  the  moments  for  large  aspect  ratios  q  can  be  obtained  at  negligible  cost 
from  a  perturbation  expansions  in  £  =  tsin0.6  Although  this  expansion  is  valid  for  all  Reynolds 
numbers,  it  disregards  the  finite  length  of  the  cylinder  and  hence  excludes  the  effect  of  resonances 
with  inertial  waves  that  depend  on  this  length. 

An  alternative  approach  to  solving  the  '  iwest-order  equations  of  such  an  expansion  has  been 
suggested  by  Hall.  Scdney,  and  Gertxt  (HSG).7  This  approach  expands  velocity  components  and 
pressure  in  a  series  of  products  of  trigonometric  functions  in  axial  direction  and  radial  “eigenfunc- 
uons”  that  satisfy  homogeneous  boundary  conditions  at  the  side  wall.  The  coefficients  of  the  series 
can  be  found  from  the  boundary  conditions  at  the  end  wall  by  collocation  or  least  squares  methods. 
Since  only  the  lineanzal  equations  are  solved,  the  HSG  method  provides  yaw  and  roll  moments  but 
the  pitch  moment  requires  solving  for  nonlinear  terms.  While  the  nonlinear  extension  may  be  a 
matter  of  time,  the  shortcoming  of  the  method  is  in  the  numerical  determination  of  eigenvalues  and 
eigenfunctions.  Practical  application  is  restricted  to  the  range  up  to  Re  =  1000  with  CPU  times  of  10 
-  1800  seconds  per  solution  (VAX  8600). 

For  given  parameters,  the  eigenvalues  for  the  HSG  expansion  are  obtained  by  iterative  solution 
of  a  sixth-order  complex  system  of  ordinary  differential  equations  and  arc  difficult  to  find.  Good  ini¬ 
tial  guesses  arc  required  for  the  iteration  to  converge.  This  problem  is  currently  overcome  by  precal¬ 
culating  voluminous  tables  for  interpolation  of  the  initial  estimates.  The  generation  of  these  tables 
requires  approximately  40  hours  CPU  time  on  a  Cray  (Murphy  1988,  personal  communication). 

In  the  following,  wc  discuss  an  alternative  to  the  HSG  method.  We  first  outline  what  quantities 
are  needed  for  the  calculation  of  the  moments.  Wc  derive  a  single  sixth-order  partial  differential 
equation  for  the  axial  velocity  component.  This  equation  can  be  solved  by  using  axial  or  radial 
eigenfunctions.  The  eigen  functions  arc  given  in  closed  form  and  the  eigenvalues  arc  determined  by 
numerically  solving  a  closed-form  characteristic  equation.  The  partial  differential  equation  in  two 
variables  also  juovides  a  new  basis  for  efficient  solution  by  spectral  method-.. 

2.  Calculation  of  the  Moments 

For  the  moments,  wc  use  cartesian  coordinates  x.y,  z,  where  z  is  the  spin  (or  cylinder)  axis 
wluie  t  is  normal  to  z  and  coplanar  with  spin  axis  and  nutation  axis.  Velocity  and  pressure  fields, 
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however,  are  more  conveniently  expressed  in  cylindrical  coordinates  r,  z.  The  nondimensional 
moments  can  be  obtained  from  the  following  relations:8 

ti  2n  1 

Mz  =  2~cos9  J  |  |  vi  r2cos $dr  d$dz  ,  (1) 

M,  =  Mx  tan0  ,  (2) 

n  2«  i  n  2*  i 

-  2— -cosB  j  |  |  vj  r2sin$dr  dfydz  +  — sin0  J  j  |  v^r2dr  dQdz  .  (3) 

The  reference  moment  is  pco2a5  and  v  =  (v,,  v^,  vt)  is  the  velocity  vector  of  the  deviation  from 
rigid-body  motion.  We  represent  the  velocity  field  by  the  Fourier  series 

v(r,  4>.  z)  =  £  V„(r,  z)  e"1*  ,  v„  =  (u„,  v„,  wn)  .  (4) 

After  performing  the  integrations  over  4>  in  eqs.  (1)  and  (3).  it  is  obvious  that  the  calculation  of  the 

moments  requires  only  the  two  Fourier  components  w,(r,  z)  and  v0(r,z).  From  perturbation 
expansion  in  e  =  xsinB  6  it  becomes  immediately  clear  that  w,  is  of  order  0(e)  and  hence  the  dom¬ 
inating  component  directly  related  to  the  periodic  forcing  with  xr  =  -  ecos<i>.  In  contrast,  v0  is  an 
aperiodic  (streaming)  component  and  hence  of  order  0(e2).  Since  e  is  small  in  cases  of  practical 
interest,  we  concentrate  in  the  following  on  determining  w ,  from  the  linearized  equations,  leaving  the 
nonlinear  corrections  for  future  investigation. 


3.  Linearized  Equations  for  the  Axial  Velocity 

We  write  the  Navier-Stokes  equations  for  the  deviation  v,  p  from  rigid  body  motion  in  the 
nutating  (aeroballistic)  coordinate  system.9  Linearization  in  e  =  (£2/o>)sin0  provides 

V  •  v  =  0  ,  (5) 

v  +  2t  x  v  +  Vp - —  72v  =  -  2 rxre*  ,  (6) 

Re 

where 


^  =  -3-tJLA+  '  JL  +  JL 

dr 2  r  dr  r 2  d<j)2  dz~ 

and 


(7) 


t  =  (' xr ,  1  +  xt )  .  x.  - -sinGcos-f)  ,  x6  =  — —  sinQ  sind>  ,  x.  =  — cosO  .  (8) 

to  to  to 

V.'c  take  the  curl  of  the  momentum  equations  (6)  to  obtain  the  vonicity  equations  that  arc  free  of  the 
pressure.  We  further  take  the  curl  of  the  vonicity  equations  and  apply  d/dd  -  {\/Re)'V2  to  the  result¬ 
ing  equations.  Since  x,  and  x0  arc  functions  of  <5  only,  wc  obtain  the  following  homogeneous  system 
of  partial  differential  equations  for  v: 


Re 


1  t  v6y  -  4-v‘«  +  v:*  +  4(i  +  x,  )!t4  =  o . 


Re  0<t> 


dz‘ 


(9) 


Wc  introduce  the  Fourier  scries  (4)  for  v  and  the  scaled  variables 
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»v 


I  =  EW  ,  f  =  qr  ,  z  =  qz  ,  q  =  { 1  +  i )'fRe/2  , 


(10) 


and  obtain  for  w 


where 


W  -  2^4w  +  V2w  -  4(1  +  x.  )2——r  =  0  , 

322 


V- jL  +  ia__l  +  JL 

"  ‘2  322 


3r' 


3r 


(11) 


(12) 


Equation  (11)  is  a  single  sixth-order  partial  differential  equation  for  w(r,  2).  The  derivation  of  the 
boundary  conditions  at  the  end  walls  and  side  wall  is  tedious  and  will  be  given  elsewhere.  The  boun¬ 
dary  conditions  at  the  end  walls  2  =  ±  <?T|  are 


w  =  0  ,  — =0,  ~  +  (2V,2  -  1)-^  =  0  , 
32  3 2s  32 3 


with 


v2  =  JL  +  l-i  -  -L  . 

3 r2  f  3 r  r 2 

Tne  boundary  conditions  at  the  side  wall  r  -  q  take  the  form 

u»  =  0  , 

-i(v.V)  +  ^ 

3r  ^  ^2 


-  -£-V2w  - 
3r  3?  32 

3w 


3r  32  2  3?  32  4 

2  -2(1+T,)4v2w  = 


3  r?2  33w  1  rr2  2(1  +  Tt) 


2(1  +tl)(-^-^:V2w  -2 


3  r  dr 


33w 
3232 2 


♦  -Vffj 

rl  dr 


-  lv2w  +  ±V*w  +  1 V2^  =  2U-t-Ti2 
r  r  ?  3  z2  «7 


(13) 


(14) 


(15a) 


(15b) 


(15c) 


Three  of  the  boundary  conditions  (13),  (15)  contain  only  derivatives  with  respect  to  2  while  the 
other  three  involve  derivatives  in  both  variables.  In  this  torm,  the  end-wall  conditions  are  homogene¬ 
ous  while  the  axial  velocity  is  forced  by  the  side-wall  conditions  (15b,  c).  The  problem  can  also  be 
reformulated  with  inhomogeneous  end-wall  conditions  and  homogeneous  conditions  at  the  side  wall 
by  subtracting  from  w  a  multiple  of  the  closed-form  solution 


w„(r)  =  — 
3 


I  |(r  ) 
/ 1(?) 


(16) 


for  the  axial  flow  in  an  infinite  cylinder.9  / ,  denotes  the  modified  Bessel  function.  The  solution  of 
the  problem  can  be  found  by  appropriate  numerical  methods  or  by  expansions  in  spatial  eigenfunc¬ 
tions.  We  note  that  the  latter  method  belongs  to  the  class  of  spectral  methods  and  is  distinguished 
from  more  generic  methods  by  the  special  definition  of  the  expansion  functions. 


Appendix  H 


109 


4.  A  Note  on  Spectra!  Solution  with  Chebyshev  Series 

Previous  work  ^  has  shown  that  spectral  methods  with  Chebyshcv-Fourier-Chebyshev  series  in 
r,  <j>,  2.  respectively,  are  efficient  in  solving  the  Navier-Stok.es  equations  in  the  four  natural  variables 
vr,  v^,  v2 ,  and  p .  If  K .  L,  M  denote  the  degrees  of  freedom  in  the  three  space  directions,  then  the 
spectral  solution  requires  determination  of  Nns  =  4  x  K  x  L  x  M  coefficients  for  all  natural  vari¬ 
ables.  The  no-slip  boundary  conditions  are  independent  of  the  variables  and  can  be  implicitly  satisfied 
by  using  linear  combinations  of  Gicbyshcv  polynomials  for  the  velocities. 

Solving  the  governing  equations  for  the  single  function  w  with  Chebyshev-Chcbyshev  series  in 
r  and  2,  respectively,  requires  the  determination  of  only  =  (K  +  2)  x  (,Vf  +  1)  expansion 
coefficients.  The  small  increases  in  the  degrees  of  freedom  in  r  and  1  account  for  the  three  boun¬ 
dary  conditions  that  involve  derivatives  with  respect  to  both  variables;  the  other  three  conditions  can 
be  implicitly  satisfied.  The  projected  CPU  time  for  a  single  solution  is  reduced  by  the  coasiderable 
factor  (Nw!NNS)'i  =  O.Olti/L3.  These  savings  should  permit  spectral  solutions  for  Reynolds  numbers 
well  in  excess  of  Re  =  104. 

The  spectral  code  was  suggested  by  M.  Sclmi  and  is  currently  in  development.  The  high  order 
of  differentiation  demands  extreme  care  for  preventing  round-off  errors  in  the  implementation. 


5.  Spatial  Eigenfunctions 

For  brevity,  we  consider  here  only  the  case  of  axial  eigenfunctions  that  satisfy  the  homogeneous 
boundary  conditions  (13).  We  solve  the  governing  equation  (12)  by  separation  of  variables: 

w(r,  z )  =  F(r)G(z)  ,  (17) 

where 

V-F  =  BF  ,  F(f)  =  /,(\I>)  ,  (18) 

and  B  is  a  constant  of  separation.  Equation  (17)  with 

3 

G(:)=  £  c:  cosa.z  (19) 

i  «  l 

is  a  solution  of  equation  (12)  if  the  three  a,  are  roots  of 

d6  -  OB  -  2 )a4  +  [3fl2  -  4B  +  1  -4(1  +  x,)2]a2  -  (B3  -  2B2  +  B )  =  0  .  (20) 


The  other  three  tools  arc  redundant  by  reasons  of  symmetry  in  z .  The  boundary  conditions  at  the  end 
walls  are  sauslied  if 


(a\-a})(2B  -  1  -  d2-a2)a]tan(a,^ri)a-,tan(a;,(7r|) 

(  (u  2  a  )  OB  -  i  -  a2  -a  j2 )  a  itar.fa  ^Tifajtanfuyyr))  (21) 

+  (:t  \  )  (?.B  -  1  -  a^-a  f  )rt3tan(rj3^q)a,tar((J  j<7 r|)  -  0 

and  two  of  the  coefficients  c,  arc  properly  chosen.  Equations  (20)  and  (21)  represent  a  transcenden¬ 
tal  nonlinear  system  that  provides  an  infinite  set  of  solutions  for  the  four  unknown  quantities  B  and 
w,  and  associated  eigenfunctions  (17). 

Finding  the  eigenvalue  quadrupels  [B .  a  a~,  a-s)  is  a  nontrivial  task,  and  details  will  be 
descried  elsewhere.  The  relatively  simple  form  of  cqs.  (20)  and  (21),  however,  permits  application 


Append i x  h 


110 


of  analytical  tools  unavailable  in  the  numerical  approach.7  Equation  (20)  has  four  isolated  roots  with 
real  B .  Only  the  double  root  B  =  1  (with  a,  =  0)  satisfies  eq.  (21)  and  provides  an  isolated  eigen¬ 
function.  For  sufficiently  large  Reynolds  numbers  (and  q),  it  can  be  shown  that  one  of  the  a, ’s, 
say.  must  be  located  near  the  diagonal  arg{a s)  =  -  tc/4.  We  observe  that  eq.  (20)  can  be  interpreted 
as  determining  cither  the  a,  for  given  B  or  three  solutions  5;  ,  j  =  1,2,3  for  any  given  a.  Varying 
a  along  the  diagonal  -  it/4  provides  three  branches  as  the  locus  of  the  Bj,  two  of  them,  B  t  and  B  2, 
originating  at  B  =  1,  the  third,  53.  at  B  =0,  as  shown  in  figure  1.  A  small  correction  to  these 
branches  and  the  position  of  the  eigenvalues  along  the  branches  is  provided  by  eq.  (21).  The  laiger 
lire  Reynolds  number,  the  deriscr  is  the  spacing  of  the  eigenvalues. 

Given  a  triplet  £; ,  we  obtain  for  each  £,  the  associated  values  of  a;t  ,  i  =  1,2,3,  where  the 
a;1  arc  located  in  a  small  neighborhood  near  the  diagonal  -  rc/4.  Figure  2  shows  the  curves  gen¬ 
erated  by  the  first  50  eigenvalues  for  Re  =  20,  T|  =  4.368,  x  =  0.16667,  and  9  =  20.  Full,  dashed, 
and  dotted  lines  indicate  the  three  families  associated  with  the  Bj  in  figure  1.  At  higher  Reynolds 
numbers,  the  lines  with  arg  (<z;1)  =  -  rt/4  are  hardly  distinguishable.  Since  each  quadrupel 
(Bj  \  u;1,  a;3)  defines  one  eigenfunction,  these  functions  arc  ordered  into  triplets  with  approxi¬ 

mately  the  same  ,.  The  triplets  can  be  ordered  by  increasing  lan  I.  We  finally  obtain  the  ordered 
eigenvalue  quadrupels  (£„,;  aHjX,  a„/2.  anji),  n  =  1.2,  •  •  •<»,  with  la„nl  <  1  ^<„+i)i  1 1  - 

The  complex  eigenfunctions  arc  oscillatory,  and  the  number  of  zeroes  of  real  and  imaginary  part 
increases  with  r. .  Analysis  of  this  oscillatory  behavior  suggests  accounting  for  the  isolated  eigen¬ 
function  with  5=1.  At  larger  Reynolds  numbers,  the  eigenfunctions  exhibit  a  pronounced 
boundary-layer  character  near  z  =  rp 

6.  Numerical  Solution  for  the  Axial  Velocity 

With  given  eigenvalues  and  eigenfunctions,  the  axial  velocity  can  be  represented  in  the  form 

W(f.z)=  £  £  AvFv(r)Gv(z)  ,  (22) 

/i  =  l  j  =•  l 

where  the  yet  unknown  constants  arc  determined  by  the  boundary  conditions  (15)  at  the  side  wall, 
r  =  q.  As  Hall,  Scdney,  and  Gerber,7  we  find  these  constants  numerically  by  using  a  collocation 
method  or  the  method  of  least  squares  and  calculate  the  associated  yaw  moment  Ml.  The  conver¬ 
gence  of  the  expansion  senes  in  terms  of  Mz  is  shown  in  figure  3.  A  relatively  large  number  of 
eigenfunctions  is  necessary  to  obtain  the  accurate  result  within  an  error  of  a  few  percent.  The  accu¬ 
racy  improves  rapidly  and  systematically  as  more  eigenfunctions  arc  taken  into  account. 

The  method  is  currently  operational  up  to  Re  =  2000.  We  have  calculated  the  eigenvalues  for 
Reynolds  numbers  as  high  as  Re  =  10A.  The  calculation  of  the  eigenfunctions  and  expansion 
cool .’ -icnis,  however,  still  suffers  from  problems  with  large  numbers  caused  by  the  exponential 
behavior  of  tngonometnc  and  Bessel  functions  of  large  complex  arguments  Woik  is  conducted  to 
circumvent  these  problems  by  proper  scaling. 

7.  Summary 

We  have  developed  a  new  analytical  formulation  to  obtain  the  liquid  moments  in  the  linear 
approximation  by  solving  a  single  partial  differential  equation  for  the  axial  velocity.  This  equation 
can  lx-  solved  by  standard  spectral  methods  or  by  use  of  expansions  in  spatial  eigenfunctions  In 
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contrast  to  earlier  work,7  we  obtain  the  eigenfunctions  in  closed  form,  and  the  eigenvalues  are 
governed  by  a  closed-form  characteristic  equation.  Our  formulation  reveals  the  structure  of  the 
eigenvalue  spectrum  and  permits  direct  calculation  of  the  eigenvalues  without  costly  lookup  tables. 
The  computer  time  for  obtaining  the  spatial  eigenfunctions  is  of  the  order  of  seconds  and  hence  negli¬ 
gible.  The  calculation  of  the  expansion  coefficients  requires  essentially  the  time  for  solving  an  alge¬ 
braic  system  with  3 N  complex  unknowns,  where  N  depends  on  the  desired  accuracy  and  increases 
with  the  Reynolds  number.  The  computational  efficiency  of  the  spatial  eigenfunction  expansion  is 
comparable  to  the  spectral  Navier-Stokcs  solver  of  Herbert  and  Li.*-  We  expect  higher  efficiency 
from  solving  for  the  axial  velocity  by  standard  spectral  methods. 
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Figure  1.  The  three  branches  of  B(a)  according 
to  eq.  (20)  for  arg(a)  =  -rc/4  generated  by  the 
first  150  eigenvalues  for  Re  =  20,  t|  ■=  4.368, 
x  -  0.16667,  and  0  =  203. 


Figure  2.  The  tlirce  families  of  three  eigenvalues 
a;(,  i  =  1,2,3,  associated  with  the  5y  of  figure 
1.  Full,  dashed,  and  dotted  line  t-  ;  for 
j  =  1,2,3,  respectively. 
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Figure  3.  Convergence  of  die  values  of  Mx  obtained  by  collocation  method  (full 
line)  and  metnod  of  least  squares  (dotted  line)  with  the  number  3 N  of  eigenfunctions. 
Re  =  20,  11  -  4.368,  x  =  0.16667,  and  0  =  20°. 


Appendix  d 


113 


